Added more models. The previous commit did not work

Author: fabriciocravo

.gitignore

@@ -15,7 +15,6 @@ text_script.py
 minority_consensus_event_figures
 linkedin_model.py
 isolate_test.py
-docs
 mobspy_jounal_models
 parameter_problem.py
 plot_config_donor.json
@@ -26,4 +25,5 @@ testing_tight_layout.png
 __pycache__/
 test_plot_images/*
 !test_plot_images/read_me.txt
-!docs/example_models/jounal_models/*
+docs/*
+!docs/example_models/

docs/example_models/journal_models/BioCRNpyler_1.py

@@ -0,0 +1,10 @@
+from mobspy import *
+
+A, B, C, D = BaseSpecies()
+
+A >> 2*B [3]
+B >> C + D [1.4]
+
+S = Simulation(A | B | C | D)
+print(S.compile())
+

docs/example_models/journal_models/BioCRNpyler_2.py

@@ -0,0 +1,51 @@
+from mobspy import *
+
+Promoter, Start_Positions, Tet, Mortal = BaseSpecies()
+Promoter.inactive, Promoter.active
+Ribo, RNA_Poly = New(Start_Positions)
+P2, Ptet = New(Promoter)
+Mrna_P2, Mrna_Ptet, GFP, RFP, CFP, GFP_F_RFP = New(Mortal)
+
+# Death here does not consider a compound for degradation.
+Mortal >> Zero [1]
+
+# Promoter activation - only Ptet is activated. P2 is always active
+Rev[Ptet.inactive + Tet >> Ptet.active][1, 1]
+
+# Read expresses the reading of RNA and DNA - both transcription and translation
+# For translation the RNA Polymerase is the reader, for protein it is the Ribossome
+def Read(Pro, R, strand):
+    sp = 'started_' + strand[0][0] # sp stands for start position
+    Start_Positions.c(sp)
+    rate = [lambda r1, r2: 2 if r1.active else 1, 1]
+    # From free position to bounded to a site
+    Rev[Pro + R.c('free_' + str(R)) >> Pro + R.c(sp).c('at_' + strand[0][0])][rate]
+    next_location = strand[1:]
+    # Movement of the reader
+    for (location, Product), (next_l, _) in zip(strand, next_location):
+        R.c(sp).c('at_' + location) >> R.c(sp).c('at_' + next_l) + Product [1]
+    # Remove reader from final location
+    R.c(sp).c('at_' + next_location[-1][0]) >> R.c(sp).c('free_' + str(R)) [1]
+
+# Zero produces nothing
+Read(Ptet, RNA_Poly, [('ptet_dna', Mrna_Ptet), ('ptet_end', Zero)])
+Read(P2, RNA_Poly, [('p2_dna', Mrna_P2), ('p2_end', Zero)])
+Read(Mrna_Ptet, Ribo, [('gfp', GFP), ('rfp', GFP_F_RFP), ('rfp_end', Zero)])
+Read(Mrna_P2, Ribo, [('cfp',  CFP), ('cfp_end',  Zero)])
+Read(Mrna_Ptet, Ribo, [('rfp', RFP), ('rfp_end', Zero)])
+
+# P2 is set to active as it is a constitutive promoter
+model = set_counts({RNA_Poly: 100, Ribo: 100, GFP: 0, RFP: 0, CFP: 0, Ptet: 1, P2.active: 1,
+                    Tet: 100, Mrna_Ptet: 0, Mrna_P2: 0, GFP_F_RFP: 0})
+S = Simulation(model)
+print(S.compile())
+
+
+
+
+
+
+
+
+
+

docs/example_models/journal_models/BioNetGen_Comparison.py

@@ -0,0 +1,22 @@
+from mobspy import *
+
+R0, L0, A0, kon, koff, kAon, kAoff, kAp, kAdp = ModelParameters(100, 500, 100, 0.01,
+                                                                0.1, 0.01, 0.1, 0.01, 0.1)
+
+r_link = BaseSpecies()
+r_link.r_0, r_link.r_1
+R, L, A = New(r_link)
+
+R.r_0 + L.r_0 >> R.r_1 + L.r_1 [kon, koff]
+
+R.a_0 + A.r_0 >> R.a_1 + A.r_1 [kAon, kAoff]
+
+(R.a_1+ L + A.n_p).r_1 >> (R.a_1 + L + A.y_p).r_1 [lambda R: kAp*R**3]
+
+A.y_p >> A.n_p [kAdp]
+
+R(R0), L(L0), A(A0)
+S = Simulation(R, L, A)
+print(S.compile())
+
+

docs/example_models/journal_models/BioNetGen_multi_state.txt

@@ -0,0 +1,42 @@
+begin parameters
+  R0   100
+  L0   500
+  A0   100
+  kon  0.01
+  koff  0.1
+  kAon  0.01
+  kAoff  0.1
+  kAp 0.01
+  kAdp 0.1
+end parameters
+
+begin molecule types
+   R(l,a)
+   L(r)
+   A(r,Y~U~P)
+end molecule types
+
+begin seed species
+  R(l,a)   R0
+  L(r)   L0
+  A(r,Y~U)   A0
+end seed species
+
+begin reaction rules
+   R(l) + L(r) <-> R(l!1).L(r!1) kon,koff
+   R(a) + A(r) <-> R(a!1).A(r!1) kAon,kAoff
+   L().R().A(Y~U) -> L().R().A(Y~P)  kAp
+   A(Y~P) -> A(Y~U)   kAdp
+end reaction rules
+
+begin observables
+  Molecules A_P  A(Y~P)
+  Molecules A_unbound_P  A(r,Y~P)
+  Molecules A_bound_P  A(r!+,Y~P)
+  Molecules RLA_P  R().L().A(Y~P)
+end observables
+
+generate_network();
+writeSBML();
+simulate_ode({t_end=>50,n_steps=>20});
+

docs/example_models/journal_models/Bistable_Gut_Inflammation.py

@@ -0,0 +1,63 @@
+from mobspy import *
+import os
+
+# Conversion taking to much time - why?
+
+Dummy, Dilutable, Promoter, TTR, Phosphorable = BaseSpecies()
+TF, Reporter = New(Dilutable)
+TtrR, TtrS = New(Phosphorable)
+Cro, CI = New(TF)
+Pcro, PcI = New(Promoter)
+Promoter.bound, Promoter.unbound, Phosphorable.dephospho, Phosphorable.phospho
+
+# TtrS triggers the state transition, TTR indicates the presence of gut inflammation
+TtrS.dephospho + TTR >> TtrS.phospho + TTR [1/u.min]
+TtrS.phospho + TtrR.dephospho >> TtrR.phospho + TtrS.dephospho [1/u.min]
+TtrR.phospho >> TtrR.dephospho [5*1e-3*1/u.second]
+TtrR.phospho >> Cro + TtrR.phospho [50*0.02*1/u.min]
+
+# Dummy represents the introduction of a constant flow of CI in the system
+Dummy >> CI + Dummy [30*0.02*1/u.min]
+
+# Dilution and degradation
+Dilutable >> Zero [0.02/u.min]
+TF >> Zero [lambda r1: 0 if r1.is_a(CI) else 1.6e-2*1/u.min]
+
+def Expression(P, Pdt, rate_expression, rate_leaky):
+    P >> P + Pdt [lambda r1: rate_expression if r1.unbound else rate_leaky]
+
+def Repression (Prom, Rep, rate_binding, rate_unbinding):
+    Prom.unbound + 2*Rep >> Prom.bound [rate_binding]
+    Prom.bound >> Prom.unbound + 2*Rep [rate_unbinding]
+
+# Pcro produces Cro, and PcI produces CI
+Expression(Pcro, Cro, 5/u.min, 0.05/u.min), Expression(PcI, CI, 4.25/u.min, 0.05/u.min)
+# Cro represses Pcro, and cI represses PcI
+Repression(Pcro, CI, 1, 50**2), Repression(PcI, Cro, 1, 40**2)
+
+model = set_counts({Pcro: 1, PcI: 1, Cro: 0, CI: 10,
+                    Reporter: 0, TtrR: 1, TtrS: 1, TTR: 0, Dummy: 0})
+S = Simulation(model)
+S.duration = 160*u.hour
+
+# Event implementation
+with S.event_time(25*u.hour):
+    TTR(1)
+with S.event_time(60*u.hour):
+    TTR(0)
+with S.event_time(90*u.hour):
+    Dummy(1)
+with S.event_time(130*u.hour):
+    Dummy(0)
+
+S.plot_data = False
+S.unit_x, S.unit_y = u.hour, 1/u.ml
+S.step_size, S.a_tol = 0.01*u.hour, 1e-16
+S.plot_config.title, S.plot_config.title_fontsize = 'Toggle Switch', 16
+S.plot_config.figsize = (6.5, 4)
+S.plot_config.vertical_lines = [25, 60, 90, 130]
+S.plot_config.ylabel = r"Conc. (mL$^{-1}$)"
+S.plot_config.save_to = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + \
+                        '/images/Toggle_Switch/Toggle_Switch.pdf'
+S.run()
+S.plot(CI, Cro)

docs/example_models/journal_models/For_The_Trees.py

@@ -0,0 +1,41 @@
+from mobspy import *
+
+if __name__ == '__main__':
+
+    """
+        This is the Tree model from the paper
+        We have a population of Trees
+        The Trees can die, age, have different colors and be in two different forests
+        The colors can change randomly from time to time
+        All old Trees can reproduce, but the Three is born green and young 
+    """
+    Age, Mortal, Colored, Location = BaseSpecies()
+    Colored.green, Colored.yellow, Colored.brown
+    Location.dense, Location.sparse
+    Age.young >> Age.old[1 / 10 / u.year]
+    Mortal >> Zero[lambda r1: 0.1 / u.year if r1.old else 0]
+    Tree = Age * Colored * Mortal * Location
+
+    # replication
+    Tree.old >> Tree + Tree.young[0.1 / u.year]
+
+    # competition
+    Tree.dense.old + Tree.dense.young >> Tree.dense.old[1e-10 * u.decimeter ** 2 / u.year]
+
+    # reproduction
+    bf = 1e-10 * u.decimeter ** 2 / u.year
+    rep_r = lambda t1, t2: 5 * bf if (Location(t1) == Location(t2) and Colored(t1) == Colored(t2)) else bf
+    2 * Tree >> 2 * Tree + Tree.yound[rep_r]
+
+    # color cycling
+    colors = ['green', 'yellow', 'brown']
+    for color, next_color in zip(colors, colors[1:] + colors[:1]):
+        Tree.c(color) >> Tree.c(next_color)[10 / u.year]
+
+    # initial conditions
+    Tree.dense(50), Tree.dense.old(50), Tree.sparse(50), Tree.sparse.old(50)
+    MySim = Simulation(Tree)
+    MySim.run(volume=1 * u.meter ** 2, unit_x='year',
+              duration=100 * u.years, repetitions=3, output_concentration=False,
+              simulation_method='stochastic', save_data=False, plot_data=False)
+    MySim.plot_stochastic(Tree.dense, Tree.sparse, Tree.brown)

docs/example_models/journal_models/Kappa_comp_v1.py

@@ -0,0 +1,34 @@
+from mobspy import *
+
+# We start with all base species
+Link1, Link2, Phos_2 = BaseSpecies()
+Link1.nl_1, Link1.l_1, Link2.nl_2, Link2.l_2
+Phos_2.not_phos, Phos_2.phos_1, Phos_2.phos_2
+
+# A and B can be linked
+A = Link1*Link2
+B = New(Link1)
+# C can be phosporolized twice
+C = Phos_2*Link2
+
+# Rev stands for reversible reaction
+Rev[A.nl_1 + B.nl_1 >> A.l_1 + B.l_1][1e-4, 0.1]
+
+# A binds to C
+A.nl_2.l_1 + C.nl_2.not_phos >> A.l_2.l_1 + C.l_2.not_phos [1e-4]
+
+# C is phosphoralized releases A
+C.l_2.not_phos + A.l_2.l_1 >> C.nl_2.phos_1 + A.nl_2.l_1 [1]
+
+# C phosporalized binds to unbound A
+A.l_2.nl_1 + C.nl_2.phos_1 >> A.l_2.nl_1 + C.l_2.phos_1 [1e-4]
+
+# Final C site is modified
+A.l_2.nl_1 + C.l_2.phos_1 >> A.nl_2.nl_1 + C.nl_2.phos_2 [1]
+
+# initial conditions
+A(1000), B(1000), C(10000)
+S = Simulation(A | B | C)
+S.compile()
+
+

docs/example_models/journal_models/Logical_XOR_Gate.py

@@ -0,0 +1,89 @@
+from mobspy import *
+import seaborn
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+max_plas, max_pbad, max_ptet  = (1, 1, 1)
+
+Mortal, Location = BaseSpecies()
+Location.c1, Location.c2, Location.c3, Location.c4
+PBad, PTet, Pcl, PLas, Movable = New(Location)
+Ara, aTc, Cl, YFP, AHL, OC12 = New(Mortal*Movable)
+n = 50
+# n = 10 - for testing
+# Parametric sweeps for XOR gate analysis
+ara_entrace_rate, atc_entrace_rate = ModelParameters([float(x) for x in np.linspace(0, 35, n)],
+                                                     [float(x) for x in np.linspace(0, 2000, n)])
+
+# Death and movement reactions
+Mortal >> Zero[1]
+for x, y in zip(['c1', 'c1', 'c2', 'c3'], ['c2', 'c3', 'c4', 'c4']):
+    Movable.c(x) >> Movable.c(y) [0.1]
+
+# Represents the entrance of aTc and Ara in the system
+def diffusion_in_cell(Molecule, rate, locations):
+    for l in locations:
+        Zero >> Molecule.c(l) [rate]
+diffusion_in_cell(Ara, ara_entrace_rate, ['c1', 'c2'])
+diffusion_in_cell(aTc, atc_entrace_rate, ['c1', 'c3'])
+
+# Promoter activation function
+def promoter_activation(P, Ligand, Protein, tf_max, n,
+                        K_d, protein_production_rate, locations):
+    tf_linked = lambda lig: tf_max * lig ** n / (lig ** n + K_d ** n)
+    tf_free = lambda lig: tf_max - tf_max * lig ** n / (lig ** n + K_d ** n)
+    pr = lambda r1, r2: protein_production_rate(r1, tf_linked(r2), tf_free(r2))
+    for l in locations:
+        with Location.c(l):
+            P + Ligand >> P + Ligand + Protein [pr]
+
+# Inverter for each nor gate
+def inverter_wire(P, R, Signal, locations):
+    rate_f = lambda r1, r2: 181*r1*350/(1 + 350 + 15*r2 + 50*r2 + 15*50*0.18*r2**2)
+    for l in locations:
+        with Location.c(l):
+            P + R >> P + R + Signal[rate_f]
+
+# Custom buffer for clear visibility
+def buffer(L, Signal, l, n, K):
+    with Location.c(l):
+        L >> L + Signal [lambda r: 30*r**n/(r**n + K**n)]
+
+# Each promoter expression is written here and assign to the promoter function
+pbad_p_rate = lambda r1, r2, r3: 765*r1*(0.009 + 37.5*r2)/(1 + 0.009 + 37.5*r2 + 3.4*r3)
+promoter_activation(PBad, Ara, Cl, max_pbad, 2.8, 90, pbad_p_rate, ['c1', 'c2'])
+ptet_p_rate = lambda r1, r2, r3: 300*r1*350/(1 + 350 + 2*160*r3 + 160**2*r3**2)
+promoter_activation(PTet, aTc, Cl, max_ptet, 1.0, 250, ptet_p_rate, ['c1', 'c3'])
+plas_p_rate = lambda r1, r2, r3: 69*r1*(0.002 + 100*r2)/(1 + 0.002 + 100*r2)
+promoter_activation(PLas, AHL, Cl, max_plas, 1.4, 0.2, plas_p_rate, ['c2', 'c3'])
+# Inverter wires in all the gates
+inverter_wire(Pcl, Cl, AHL, ['c1']), inverter_wire(Pcl, Cl, OC12, ['c2', 'c3'])
+# Buffer for YFP
+buffer(OC12, YFP, 'c4', 4, 0.04)
+
+# Count assignment and Sim
+model = set_counts({Ara: 0, aTc: 0, Cl: 0, YFP: 0, PBad.c1: 1, PBad.c2: 1, PTet.c1: 1, PTet.c3: 1,
+                    Pcl.c1: 1, Pcl.c2: 1, Pcl.c3: 1, PLas.c2: 1, PLas.c3: 1, AHL: 0, OC12: 0})
+S = Simulation(model)
+S.duration = 100
+S.plot_data = False
+S.run()
+
+matrix = []
+line = []
+for i, r in enumerate(S.results[YFP]):
+    line.append(r[-1])
+    if (i + 1) % n == 0:
+        matrix.append(line)
+        line = []
+
+ax = seaborn.heatmap(matrix)
+
+# Invert the y-axis
+ax.invert_yaxis()
+
+plt.title('XOR Gate Heatmap', fontsize=18)
+plt.xlabel('Ara (a.u.)', fontsize=14)
+plt.ylabel('aTc (a.u.)', fontsize=14)
+plt.show()