diff --git a/Tests/Unit/test_chemicalReactionNetwork.py b/Tests/Unit/test_chemicalReactionNetwork.py index 6ba2e114..8d5256b6 100644 --- a/Tests/Unit/test_chemicalReactionNetwork.py +++ b/Tests/Unit/test_chemicalReactionNetwork.py @@ -228,7 +228,7 @@ def test_initial_condition_vector(self): def test_get_all_species_containing(self): # test that the species arument must be Species object with self.assertRaisesRegex( - ValueError, 'species argument must be an instance of Species!' + ValueError, 'species argument must be an instance of Species' ): self.crn.get_all_species_containing(species=self.species_list) diff --git a/Tests/Unit/test_component_basic.py b/Tests/Unit/test_component_basic.py index 55caf0d5..f02de3e0 100644 --- a/Tests/Unit/test_component_basic.py +++ b/Tests/Unit/test_component_basic.py @@ -11,7 +11,11 @@ def test_enzyme(): products = 'P1' enzyme = 'E1' - e = Enzyme(enzyme=enzyme, substrates=substrates, products=products) + e = Enzyme( + enzyme=enzyme, substrates=substrates, products=products, + ) + e.default_mechanism=None # override failsafe + assert any([s.name == 'S1' for s in e.substrates]) assert any([p.name == 'P1' for p in e.products]) assert enzyme == e.enzyme.name diff --git a/Tests/Unit/test_component_membrane.py b/Tests/Unit/test_component_membrane.py index 4a2f9126..99db5a29 100644 --- a/Tests/Unit/test_component_membrane.py +++ b/Tests/Unit/test_component_membrane.py @@ -21,6 +21,7 @@ def test_DiffusibleMolecule(): assert dm.get_species().name == 'DP' + dm.default_mechanism = None # override failsafe with pytest.raises( KeyError, match='Unable to find mechanism of type diffusion in Component', @@ -46,6 +47,7 @@ def test_IntegralMembraneProtein(): assert imp.get_species().name == 'MP1' + imp.default_mechanism = None # override failsafe with pytest.raises( KeyError, match='Unable to find mechanism of type membrane_insertion in Component', @@ -70,6 +72,7 @@ def test_MembraneChannel(): assert mc.get_species().name == 'IMP1' + mc.default_mechanism = None # override failsafe with pytest.raises( KeyError, match='Unable to find mechanism of type transport in Component', @@ -94,6 +97,7 @@ def test_MembranePump(): assert mp.get_species().name == 'MPump1' + mp.default_mechanism = None # override failsafe with pytest.raises( KeyError, match='Unable to find mechanism of type transport in Component', @@ -125,6 +129,7 @@ def test_MembraneSensor(): assert ms.get_species().name == 'MSensor1' + ms.default_mechanism = None # override failsafe with pytest.raises( KeyError, match='Unable to find mechanism of type membrane_sensor in Component', diff --git a/biocrnpyler/__init__.py b/biocrnpyler/__init__.py index 1c443472..16504f57 100644 --- a/biocrnpyler/__init__.py +++ b/biocrnpyler/__init__.py @@ -11,6 +11,12 @@ # All utilities from .utils import * +from .utils.plotting import ( + plot_all_species_containing as plot_all_species_containing, +) +from .utils.plotting import ( + plot_gene_expression_data as plot_gene_expression_data, +) try: from ._version import version as __version__ diff --git a/biocrnpyler/components/basic.py b/biocrnpyler/components/basic.py index 30088056..014533f8 100644 --- a/biocrnpyler/components/basic.py +++ b/biocrnpyler/components/basic.py @@ -6,6 +6,9 @@ from ..core.component import Component from ..core.reaction import Reaction from ..core.species import Complex, Species +from ..mechanisms.binding import One_Step_Binding +from ..mechanisms.enzyme import MichaelisMenten +from ..mechanisms.metabolite import OneStepPathway class DNA(Component): @@ -402,7 +405,17 @@ def __init__( else: self.products.append(self.set_species(p)) - Component.__init__(self=self, name=name, **kwargs) + if precursors is not None or products is not None: + default_mechanism = OneStepPathway() + else: + default_mechanism = None + + Component.__init__( + self=self, + name=name, + default_mechanism=default_mechanism, + **kwargs, + ) def get_species(self) -> Species: """Get the metabolite species. @@ -602,7 +615,12 @@ def __init__( if name is None: name = self.species.name - Component.__init__(self=self, name=name, **kwargs) + Component.__init__( + self=self, + name=name, + default_mechanism=One_Step_Binding(), + **kwargs, + ) def get_species(self) -> List[Species]: """Get the complex species. @@ -773,7 +791,12 @@ def __init__( self.substrates = substrates self.products = products - Component.__init__(self=self, name=self.enzyme.name, **kwargs) + Component.__init__( + self=self, + name=self.enzyme.name, + default_mechanism=MichaelisMenten(), + **kwargs, + ) @property def substrates(self) -> List: diff --git a/biocrnpyler/components/combinatorial_conformation.py b/biocrnpyler/components/combinatorial_conformation.py index ed49a49f..df15dc55 100644 --- a/biocrnpyler/components/combinatorial_conformation.py +++ b/biocrnpyler/components/combinatorial_conformation.py @@ -140,8 +140,12 @@ def __init__( if name is None: name = str(self.internal_polymer.name) - Component.__init__(self, name=name, **kwargs) - self.add_mechanism(One_Step_Reversible_Conformation_Change()) + Component.__init__( + self, + name=name, + default_mechanism=One_Step_Reversible_Conformation_Change(), + **kwargs, + ) # Helper function to assert the correct class type def _assert_conformation(self, states, input_name='states'): diff --git a/biocrnpyler/components/membrane.py b/biocrnpyler/components/membrane.py index 026423df..61430a8a 100644 --- a/biocrnpyler/components/membrane.py +++ b/biocrnpyler/components/membrane.py @@ -7,6 +7,14 @@ from ..core.compartment import Compartment from ..core.component import Component from ..core.species import Species +from ..mechanisms.signaling import Membrane_Signaling_Pathway_MM +from ..mechanisms.transport import ( + Facilitated_Transport_MM, + Membrane_Protein_Integration, + Primary_Active_Transport_MM, + Simple_Diffusion, + Simple_Transport, +) class DiffusibleMolecule(Component): @@ -108,7 +116,11 @@ def __init__( name = self.substrate.name + '_' + self.substrate.compartment.name Component.__init__( - self=self, name=name, attributes=attributes, **kwargs + self=self, + name=name, + attributes=attributes, + default_mechanism=Simple_Diffusion(), + **kwargs, ) def get_species(self): @@ -334,7 +346,12 @@ def __init__( + self.membrane_protein.compartment.name ) - Component.__init__(self=self, name=name, **kwargs) + Component.__init__( + self=self, + name=name, + default_mechanism=Membrane_Protein_Integration(), + **kwargs, + ) def get_species(self): """Get the membrane protein species before insertion. @@ -504,7 +521,7 @@ def __init__( integral_membrane_protein = self.set_species( integral_membrane_protein, material_type='protein', - attributes='Passive' if direction is None else direction, + attributes=['Passive'] if direction is None else direction, ) self.integral_membrane_protein = integral_membrane_protein @@ -568,7 +585,12 @@ def __init__( + self.integral_membrane_protein.compartment.name ) - Component.__init__(self=self, name=name, **kwargs) + Component.__init__( + self=self, + name=name, + default_mechanism=Simple_Transport(), + **kwargs, + ) def get_species(self): """Get the integral membrane protein species. @@ -782,14 +804,14 @@ def __init__( self.membrane_pump = self.set_species( membrane_pump, material_type='protein', - attributes='Passive', + attributes=['Passive'], ) self.membrane_pump.ATP = ATP else: self.membrane_pump = self.set_species( membrane_pump, material_type='protein', - attributes=direction, + attributes=[direction], ) self.membrane_pump.ATP = ATP if direction == 'Importer': @@ -812,19 +834,19 @@ def __init__( self.integral_membrane_protein = self.set_species( membrane_pump, material_type='protein', - attributes='Passive', + attributes=['Passive'], ) elif membrane_pump.attributes[0] == 'Exporter': self.membrane_pump = self.set_species( membrane_pump, material_type='protein', - attributes='Exporter', + attributes=['Exporter'], ) elif membrane_pump.attributes[0] == 'Importer': self.membrane_pump = self.set_species( membrane_pump, material_type='protein', - attributes='Importer', + attributes=['Importer'], ) self.energy = self.set_species( 'ATP', @@ -867,7 +889,18 @@ def __init__( + self.membrane_pump.compartment.name ) - Component.__init__(self=self, name=name, **kwargs) + # Determine the default mechanism + if 'Passive' in self.membrane_pump.attributes: + default_mechanism = Primary_Active_Transport_MM() + else: + default_mechanism = Facilitated_Transport_MM() + + Component.__init__( + self=self, + name=name, + default_mechanism=default_mechanism, + **kwargs, + ) def get_species(self): """Get the membrane pump protein species. @@ -1137,7 +1170,12 @@ def __init__( + self.membrane_sensor_protein.compartment.name ) - Component.__init__(self=self, name=name, **kwargs) + Component.__init__( + self=self, + name=name, + default_mechanism=Membrane_Signaling_Pathway_MM(), + **kwargs, + ) def get_species(self): """Get the membrane sensor protein species. diff --git a/biocrnpyler/core/chemical_reaction_network.py b/biocrnpyler/core/chemical_reaction_network.py index 171bb008..dd711a79 100644 --- a/biocrnpyler/core/chemical_reaction_network.py +++ b/biocrnpyler/core/chemical_reaction_network.py @@ -604,7 +604,7 @@ def initial_condition_vector( return x0 def get_all_species_containing( - self, species: Species, return_as_strings=False + self, species: Union[Species, str], return_as_strings=False ): """Find all species (complexes) that contain a given species. @@ -613,8 +613,9 @@ def get_all_species_containing( Parameters ---------- - species : Species - The species to search for within other species. + species : Species or str + The species to search for within other species. If given as a + string, returns all species with that string as part of the name. return_as_strings : bool, default=False If True, returns species as string representations. If False, returns actual Species objects. @@ -645,13 +646,15 @@ def get_all_species_containing( """ return_list = [] - if not isinstance(species, Species): + if not isinstance(species, (Species, str)): raise ValueError( - "species argument must be an instance of Species!" + "species argument must be an instance of Species or str." ) for s in self._species: - if species in s.get_species(recursive=True): + if ( + isinstance(species, str) and species in s.name + ) or species in s.get_species(recursive=True): if return_as_strings: return_list.append(repr(s)) else: diff --git a/biocrnpyler/core/component.py b/biocrnpyler/core/component.py index 749cff5c..e264ebea 100644 --- a/biocrnpyler/core/component.py +++ b/biocrnpyler/core/component.py @@ -54,6 +54,8 @@ class Component: initial_condition_dictionary : dict, optional Dictionary mapping species (or species names) to initial concentration values for components with multiple species. + default_mechanism : Mechanism, optional + Default mechanism to use if no other mechanisms are provided. Attributes ---------- @@ -120,11 +122,15 @@ def __init__( # None, self.name):initial_concentration initial_concentration=None, initial_condition_dictionary=None, + default_mechanism=None, ): if mechanisms is None: self.mechanisms = {} + self.default_mechanism = default_mechanism else: self.mechanisms = mechanisms + self.default_mechanism = None + if isinstance(name, Species): self.name = name.name elif isinstance(name, str): @@ -484,6 +490,10 @@ def get_mechanism(self, mechanism_type, optional_mechanism=False): return self.mechanisms[mechanism_type] elif self.mixture is not None: mech = self.mixture.get_mechanism(mechanism_type) + + if mech is None and self.default_mechanism is not None: + mech = self.default_mechanism + if mech is None and not optional_mechanism: raise KeyError( f"Unable to find mechanism of type {mechanism_type} in " @@ -643,7 +653,7 @@ def get_parameter( Name of the parameter to retrieve. part_id : str, optional Part identifier for the parameter lookup key. - mechanism : str, optional + mechanism : str or Mechanism, optional Mechanism identifier for the parameter lookup key. return_numerical : bool, default=False If True, returns the numerical value. If False, returns the @@ -688,19 +698,23 @@ def get_parameter( if param is None and self.mixture is not None and check_mixture: param = self.mixture.get_parameter(mechanism, part_id, param_name) - # Finally, try the mechanism parameter database - if param is None and self.mechanisms is not None and check_mechanism: + # Finally, try mechanism parameter databases + if param is None and check_mechanism: + # Start with the component mechanism if isinstance(mechanism, Mechanism): key = mechanism.mechanism_type else: key = mechanism - if self.mechanisms.get(key, None): + if self.mechanisms is not None and self.mechanisms.get(key, None): param = self.mechanisms[key].get_parameter( mechanism, part_id, param_name ) - elif ( - self.mixture + + # Then mixture mechanisms + if ( + param is None + and self.mixture and self.mixture.mechanisms and self.mixture.mechanisms.get(key, None) ): @@ -708,6 +722,12 @@ def get_parameter( mechanism, part_id, param_name ) + # Then the passed mechanism + if param is None and isinstance(mechanism, Mechanism): + param = mechanism.get_parameter( + mechanism, part_id, param_name + ) + if param is None and not return_none: raise ValueError( "No parameters can be found that match the " diff --git a/biocrnpyler/mechanisms/binding_parameters.tsv b/biocrnpyler/mechanisms/binding_parameters.tsv index 4427719f..7bc66805 100644 --- a/biocrnpyler/mechanisms/binding_parameters.tsv +++ b/biocrnpyler/mechanisms/binding_parameters.tsv @@ -8,24 +8,41 @@ mechanism_id part_id param_name param_val units comments # # One step (and cooperative) binding # -# The first step in various binding mechanisms is the initial binding of -# the substrates. We define binding and unbinding reaction rates here -# for several different part types. These rates are used by both -# 'one_step_binding' mechanisms and ' -# Chemical complex binding -binding chemical_complex kb 100 /sec assuming 10 ms to diffuse across 1 um (characteristic cell size) -binding chemical_complex ku 10 uM/sec 90% binding +# The first step in various binding mechanisms is the initial binding +# of the substrates. We define binding and unbinding reaction rates +# here for two different situations: a chemical complex and a +# regulator (activator or repressor) binding to DNA. +# +# We use TetR and aTc as the "generic" case and use the rate constants +# from the following paper: +# +# [SK07] V. Sotiropoulos and Y. N. Kaznessis. "Synthetic +# Tetracycline-Inducible Regulatory Networks: Computer-Aided Design +# of Dynamic Phenotypes." BMC Systems Biology 1, no. 1 (2007): 7. +# https://doi.org/10.1186/1752-0509-1-7. +# +# Table 1 presents rates for a large set of reactions (using Tc instead +# of aTc) in units of (M s)^{-1} (converted to (uM s)^{-1} here by +# multiplying by 1e-6. + +# Chemical complex binding (aTc to TetR) +binding chemical_complex kb 2.00 /uM/sec SK07, Table 1, R6 +binding chemical_complex ku 0.001 /sec SK07, Table 1, R8 -# Protein-DNA binding -binding dna_protein kb 0.2 /sec TODO -binding dna_protein ku 1e-3 uM/sec TODO +# Protein-DNA binding (TetR2 to DNA) +binding dna_protein kb 2.86 /uM/sec SK07, Table 1, R23 +binding dna_protein ku 5.11e-4 /sec SK07, Table 1, R24 # Default cooperativity binding cooperativity 2 number of binding modules to bind simultaneously -# Two step cooperative binding (oligimerization) -two_step_cooperative_binding kb1 1.0 /sec forward rate constant for oligomerization (TODO) -two_step_cooperative_binding ku1 0.01 uM/sec reverse rate constant for oligomerization (TODO) -two_step_cooperative_binding kb2 1.0 /sec forward rate constant for oligomer-bindee binding (TODO) -two_step_cooperative_binding ku2 0.01 uM/sec reverse rate constant for oligomer-bindee binding (TODO) +# Two step cooperative binding: +# R1 (oligimerization): 2TetR --> TetR2) +# R2 (DNA binding): TetR2 to DNA) +two_step_cooperative_binding kb1 1.0e3 /uM/sec SK07, Table 1, Reaction 2 +two_step_cooperative_binding ku1 10 /sec SK07, Table 1, Reaction 3 +two_step_cooperative_binding chemical_complex kb2 2.00 /uM/sec SK07, Table 1, Reaction 6 +two_step_cooperative_binding chemical_complex ku2 0.001 /sec SK07, Table 1, Reaction 8 +two_step_cooperative_binding dna_protein kb2 2.86 /uM/sec SK07, Table 1, Reaction 23 +two_step_cooperative_binding dna_protein ku2 5.11e-4 /sec SK07, Table 1, Reaction 24 diff --git a/biocrnpyler/mechanisms/enzyme_parameters.tsv b/biocrnpyler/mechanisms/enzyme_parameters.tsv index 259da579..7f2647e3 100644 --- a/biocrnpyler/mechanisms/enzyme_parameters.tsv +++ b/biocrnpyler/mechanisms/enzyme_parameters.tsv @@ -6,13 +6,13 @@ mechanism_id part_id param_name param_val units comments # Basic (Michaelis-Menten) catalysis parameters -catalysis kb 100 /sec assuming 10 ms to diffuse across 1 um (characteristic cell size) -catalysis ku 10 uM/sec 90% binding -# catalysis kcat 1 +catalysis kb 100 /uM/sec assuming 10 ms to diffuse across 1 um (characteristic cell size) +catalysis ku 10 /sec 90% binding +# catalysis kcat 1 /sec # Reversible Michaelis-Menten catalysis parameters -catalysis kb1 100 /sec assuming 10 ms to diffuse across 1 um (characteristic cell size) -catalysis ku1 10 uM/sec 90% binding +catalysis kb1 100 /uM/sec assuming 10 ms to diffuse across 1 um (characteristic cell size) +catalysis ku1 10 sec 90% binding catalysis kb2 100 /sec assuming 10 ms to diffuse across 1 um (characteristic cell size) catalysis ku2 10 uM/sec 90% binding catalysis kcat_rev 0.01 diff --git a/biocrnpyler/mechanisms/txtl.py b/biocrnpyler/mechanisms/txtl.py index 7ca67309..c7db48c5 100644 --- a/biocrnpyler/mechanisms/txtl.py +++ b/biocrnpyler/mechanisms/txtl.py @@ -1983,7 +1983,7 @@ def __init__( else: raise ValueError("Fuels must be a list of Species!") - if all([isinstance(s, Species) for s in fuels]): + if all([isinstance(s, Species) for s in wastes]): self.wastes = wastes else: raise ValueError("wastes must be a list of Species!") diff --git a/biocrnpyler/mechanisms/txtl_parameters.tsv b/biocrnpyler/mechanisms/txtl_parameters.tsv index 440f9a87..2c920392 100644 --- a/biocrnpyler/mechanisms/txtl_parameters.tsv +++ b/biocrnpyler/mechanisms/txtl_parameters.tsv @@ -16,31 +16,45 @@ gene_expression kexpress 19.23 /sec transcription ktx 0.05 /sec 50 nt/s and transcript length of 1000 nt translation ktl 0.05 /sec 15 aa/s and protein length of 300 aa +# # PositiveHillTranscription and NegativeHillTranscription -transcription K 1.6e-3 uM RMM (TODO) -transcription k 3.25 /sec/uM Tx_cat from Singhal (= ktx?? TODO) +# +# We computing the various rates here based on the the following paper: +# +# [SK07] V. Sotiropoulos and Y. N. Kaznessis. "Synthetic +# Tetracycline-Inducible Regulatory Networks: Computer-Aided Design +# of Dynamic Phenotypes." BMC Systems Biology 1, no. 1 (2007): 7. +# https://doi.org/10.1186/1752-0509-1-7. +# +# Table 1 presents rates for a large set of reactions (using Tc instead +# of aTc) in units of (M s)^{-1} (converted to (uM s)^{-1} here by +# multiplying by 1e-6. Dissociation constaints are computed by dividing +# unbinding rate by the binding rate. + +transcription K 1.33e-3 uM SK01 -> BFS formulas w/ RNAP = 1.3 (PURE) +transcription k 2.5e-2 /sec hand-tuned to match MM in gfp_expression.py (TODO) transcription n 2 cooperativity -transcription kleak 5e-6 /sec RMM (TODO) +transcription kleak 0.013 /sec SK07, Table 1, R37 # Transcription_MM and Translation_MM transcription_mm ktx 0.05 /sec 50 nt/s and transcript length of 1000 nt -transcription_mm kb 4.48 /uM/sec txtlsim S3: pol_{F, lac} -transcription_mm ku 2.5E-06 /sec default_parameters.txt; source unclear (TODO) +transcription_mm kb 8.6 /uM/sec SK07, Table 1, R31 +transcription_mm ku 0.01 /sec SK07, Table 1, R32 translation_mm ktl 0.05 /sec 15 aa/s and protein length of 300 aa -translation_mm kb 0.819 /uM/sec Singhal et al. Table S2 (TODO) -translation_mm ku 0.003 /sec Singhal et al. Table S2 (TODO) +translation_mm kb 0.1 /uM/sec SK07, Table 1, R45 +translation_mm ku 100 /sec SK07, Table 1, R46 # Energy_Transcription_MM and Energy_Translation_MM energy_transcription_mm ktx 50 /sec 50 nt/s (rate will be divided by length) -energy_transcription_mm kb 4.48 /uM/sec txtlsim S3: pol_{F, lac} -energy_transcription_mm ku 2.5E-06 /sec default_parameters.txt; source unclear (TODO) +energy_transcription_mm kb 8.6 /uM/sec SK07, Table 1, R31 +energy_transcription_mm ku 0.01 /sec SK07, Table 1, R32 energy_transcription_mm kb_ntps 10 /uM/sec fast binding of NTPs to RNAP:DNA complex (TODO) energy_transcription_mm ku_ntps 1e-6 /sec slow unbinding of NTPs to RNAP:DNA complex (TODO) energy_transcription_mm length 1000 nt default length of gene (for resource usage) energy_translation_mm ktl 15 /sec 15 aa/s (rate will be divided by length) -energy_translation_mm kb 0.819 /uM/sec Singhal et al. Table S2 -energy_translation_mm ku 0.003 /sec Singhal et al. Table S2 +energy_translation_mm kb 0.1 /uM/sec SK07, Table 1, R45 +energy_translation_mm ku 100 /sec SK07, Table 1, R46 energy_translation_mm kb_fuel 10 /uM/sec fast binding of NTPs to RNAP:DNA complex (TODO) energy_translation_mm ku_fuel 1e-6 /sec slow unbinding of NTPs to RNAP:DNA complex (TODO) energy_translation_mm length 300 aa default length of a protein (for resource uasge) diff --git a/biocrnpyler/mixtures/pure.py b/biocrnpyler/mixtures/pure.py index 8a549f49..dce886bb 100644 --- a/biocrnpyler/mixtures/pure.py +++ b/biocrnpyler/mixtures/pure.py @@ -3,9 +3,246 @@ from ..components.basic import Metabolite, Protein from ..core.mixture import Mixture -from ..mechanisms.binding import One_Step_Binding -from ..mechanisms.enzyme import MichaelisMenten -from ..mechanisms.txtl import Energy_Transcription_MM, Energy_Translation_MM +from ..mechanisms.txtl import ( + Energy_Transcription_MM, + Energy_Translation_MM, + SimpleTranscription, + SimpleTranslation, + Transcription_MM, + Translation_MM, +) + + +class PURE(Mixture): + """PURE cell-free protein synthesis system with customizable mechanisms. + + A mixture that models the PURE (Protein synthesis Using Recombinant + Elements) reconstituted cell-free transcription-translation system with + customizable mechanisms. + + Parameters + ---------- + name : str, default='PURE' + Name identifier for the mixture. + include_machinery : bool, default=True + Include components and mechanisms for RNAP and ribosomes. + include_resources : bool, default=True + Include components and mechanisms for NTP and amino acid utilization. + Requires `include_machinery` to be true. + include_fuel : bool, default=True + Include fuel component and energy utilization mechanism. Requires that + `include_machinery` and `include_resources` be True. + rnap : str, default='RNAP' + Name for the RNA polymerase protein species. + ribosome : str, default='Ribo' + Name for the ribosome protein species. + ntps : str, default='NTPs' + Name for the nucleotide triphosphate species (lumped NTPs excluding + ATP). + ndps : str, default='NDPs' + Name for the nucleotide diphosphate species (lumped NDPs). + amino_acids : str, default='AAs' + Name for the amino acid species (lumped amino acids). + fuel : str, default='ATP' + Name for the primary energy carrier species (ATP). + parameter_file : str, default='mixtures/pure_parameters.tsv' + Path to file containing default parameter values for the PURE + system. + **kwargs + Additional keyword arguments passed to the parent Mixture class. + + Attributes + ---------- + rnap : Protein + RNA polymerase component. + ribosome : Protein + Ribosome component. + ntps : Metabolite + Nucleotide triphosphate metabolite component (excluding ATP). + amino_acids : Metabolite + Amino acid metabolite component. + fuel : Metabolite + Fuel metabolite component (ATP). + name : str + Name of the mixture. + + See Also + -------- + EnergyTxTlExtract : TX-TL with fuel regeneration. + TxTlExtract : TX-TL with machinery but no energy. + Energy_Transcription_MM : Mechanism for energy-consuming transcription. + Energy_Translation_MM : Mechanism for energy-consuming translation. + Mixture : Base class for all mixtures. + + Notes + ----- + This mixture automatically adds the following components: + + - RNA polymerase (RNAP) + - Ribosome + - Amino acids (lumped) + - NTPs (nucleotide triphosphates excluding ATP, lumped) + - NDPs (nucleotide diphosphates, lumped) + - Fuel (ATP for energy) + + Default mechanisms included: + + - 'transcription' : `Energy_Transcription_MM` - Michaelis-Menten + transcription with length-dependent ATP and NTP consumption + - 'translation' : `Energy_Translation_MM` - Michaelis-Menten translation + with length-dependent amino acid and ATP consumption + - 'catalysis' : `MichaelisMenten` - General Michaelis-Menten enzyme + catalysis for user-defined enzymatic reactions + - 'binding' : `One_Step_Binding` - Simple multi-species binding for + forming complexes + + Key features of this mixture: + + - Explicit modeling of PURE system components + - Length-dependent energy consumption (realistic stoichiometry) + - No fuel regeneration mechanisms (finite resource pool) + - Resource competition effects (genes compete for RNAP and ribosomes) + - Resource depletion dynamics (ATP, NTPs, amino acids deplete) + - Enzyme sequestration in complexes + - Separate tracking of ATP vs other NTPs + - Suitable for modeling batch-mode PURE reactions + + Energy model details: + + - Transcription: Consumes L NTPs and L ATPs per mRNA of length L + - Translation: Consumes L amino acids and 4L ATPs per protein of length + L (4 ATPs per amino acid reflect GTP hydrolysis during elongation) + - No regeneration: ATP, NTPs, and amino acids are consumed but not + regenerated + - Energy depletion: Expression stops when resources are exhausted + - Length parameter L: Represents gene/protein length in appropriate + units + - Lumped species: Different nucleotides lumped into NTPs, different + amino acids lumped into single species + - Separate ATP: ATP tracked separately from other NTPs for independent + energy accounting + + Differences from `EnergyTxTlExtract`: + + - No fuel regeneration pathway (no NTP regeneration from 3PGA or other + fuel sources) + - ATP modeled as separate fuel species rather than included in NTPs + - Default parameter file points to PURE-specific parameters + - Intended for modeling finite-resource batch reactions + - More realistic for in vitro PURE systems + + Common applications include: + + - PURE cell-free TX-TL systems + - Resource-limited gene expression modeling + - TX-TL system optimization with fixed resource budgets + - Batch mode TX-TL reactions + - Energy budget and resource allocation studies + - Multi-gene expression burden analysis + - In vitro synthetic biology applications + + Examples + -------- + Create a PURE mixture for GFP expression: + + >>> gfp_gene = bcp.DNAassembly( + ... name='gfp_construct', + ... promoter='pconst', + ... rbs='bcd2', + ... transcript='gfp_mrna', + ... protein='GFP' + ... ) + >>> mixture = bcp.BasicPURE( + ... name='pure_mixture', + ... components=[gfp_gene], + ... parameter_file='mixtures/pure_parameters.tsv' + ... ) + >>> crn = mixture.compile_crn() + + """ + + def __init__( + self, + include_machinery=True, + include_resources=True, + include_fuel=True, + name='PURE', + rnap='RNAP', + ribosome='Ribo', + ntps='NTPs', + ndps='NDPs', + amino_acids='AAs', + atp='ATP', + adp='ADP', + fuel='Fuel_CP', + parameter_file='mixtures/pure_parameters.tsv', + **kwargs, + ): + Mixture.__init__( + self, name=name, parameter_file=parameter_file, **kwargs + ) + + # Start with basic mechansims for transcription and translation + default_mechanisms = { + 'transcription': SimpleTranscription(), + 'translation': SimpleTranslation(), + } + + if include_fuel: + if not include_resources: + raise ValueError("include_fuel requires include_resources") + self.fuel = Metabolite(fuel) + self.adp = Metabolite(adp) + self.atp = Metabolite( + atp, precursors=[self.fuel, self.adp], products=[self.adp] + ) # fuel becomes ATP, and ATP is degraded + self.add_components([self.atp]) # includes ADP, fuel + else: + self.atp = None + self.adp = None + self.fuel = None + + if include_machinery: + self.rnap = Protein(rnap) + self.ribosome = Protein(ribosome) + self.add_components([self.rnap, self.ribosome]) + + default_mechanisms |= { + 'transcription': Transcription_MM(rnap=self.rnap.species), + 'translation': Translation_MM(ribosome=self.ribosome.species), + } + + if include_resources: + if not include_machinery: + raise ValueError( + "include_resources requires include_machinery" + ) + + # create default Components to represent cellular machinery + self.ntps = Metabolite(ntps) + self.amino_acids = Metabolite(amino_acids) + self.add_components([self.amino_acids, self.ntps]) + + # Create default TX-TL Mechanisms + default_mechanisms['transcription'] = Energy_Transcription_MM( + rnap=self.rnap.get_species(), + fuels=[self.ntps.get_species()], + wastes=[], + ) + if include_fuel: + default_mechanisms['translation'] = Energy_Translation_MM( + ribosome=self.ribosome.get_species(), + fuels=4 * [self.atp.get_species()] + + [self.amino_acids.get_species()], + wastes=4 * [self.adp.get_species()], + ) + else: + default_mechanisms['translation'] = Energy_Translation_MM( + ribosome=self.ribosome.get_species(), + fuels=[self.amino_acids.get_species()], + wastes=[], + ) + self.add_mechanisms(mechanisms=default_mechanisms, overwrite=False) class BasicPURE(Mixture): @@ -176,45 +413,12 @@ def __init__( parameter_file='mixtures/pure_parameters.tsv', **kwargs, ): - Mixture.__init__( - self, name=name, parameter_file=parameter_file, **kwargs - ) - - # create default Components to represent cellular machinery - self.rnap = Protein(rnap) - self.ribosome = Protein(ribosome) - self.ntps = Metabolite(ntps) - self.amino_acids = Metabolite(amino_acids) - self.fuel = Metabolite(fuel) - - default_components = [ - self.rnap, - self.ribosome, - self.amino_acids, - self.ntps, - self.fuel, - ] - self.add_components(default_components) - - # Create default TX-TL Mechanisms - mech_tx = Energy_Transcription_MM( - rnap=self.rnap.get_species(), - fuels=[self.ntps.get_species()], - wastes=[], - ) - mech_tl = Energy_Translation_MM( - ribosome=self.ribosome.get_species(), - fuels=4 * [self.fuel.get_species()] - + [self.amino_acids.get_species()], - wastes=[], + PURE.__init__( + self, + name=name, + parameter_file=parameter_file, + include_machinery=True, + include_resources=True, + include_fuel=True, + **kwargs, ) - mech_cat = MichaelisMenten() - mech_bind = One_Step_Binding() - - default_mechanisms = { - mech_tx.mechanism_type: mech_tx, - mech_tl.mechanism_type: mech_tl, - mech_cat.mechanism_type: mech_cat, - mech_bind.mechanism_type: mech_bind, - } - self.add_mechanisms(default_mechanisms, overwrite=None) diff --git a/biocrnpyler/mixtures/pure_parameters.tsv b/biocrnpyler/mixtures/pure_parameters.tsv index 214e86ca..553a7b2a 100644 --- a/biocrnpyler/mixtures/pure_parameters.tsv +++ b/biocrnpyler/mixtures/pure_parameters.tsv @@ -39,6 +39,12 @@ initial concentration ATP 1000 uM [ATP/GTP] - [CTP/UTP] initial concentration NTPs 4000 uM [CTP/UTP] * 4 (NTPs) initial concentration AAs 6000 uM 0.3 mM (upper bound) * 20 AAs +transcription_mm kb 4.48 /uM/sec txtlsim S3: pol_{Flac} +transcription_mm ku 2.5E-06 /sec default_parameters.txt; source unclear (TODO) + +translation_mm kb 0.819 /uM/sec Singhal et al. Table S2 +translation_mm ku 0.003 /sec Singhal et al. Table S2 + energy_transcription_mm ktx 50 /sec 50 nt/s (rate will be divided by length) energy_transcription_mm kb 4.48 /uM/sec txtlsim S3: pol_{Flac} energy_transcription_mm ku 2.5E-06 /sec default_parameters.txt; source unclear (TODO) @@ -52,3 +58,8 @@ energy_translation_mm ku 0.003 /sec Singhal et al. Table S2 energy_translation_mm kb_fuel 10 /uM/sec fast binding of NTPs to RNAP:DNA complex (TODO) energy_translation_mm ku_fuel 1e-6 /sec slow unbinding of NTPs to RNAP:DNA complex (TODO) energy_translation_mm length 300 aa default length of a protein (for resource uasge) + +# Parameters for energy regeneration +initial concentration Fuel_CP 16666 uM 50 mM small molecule * 3/10 in final mix +one_step_pathway ATP_production k 0.02 /uM/sec alpha_atp, Singhal et al. Table S2 (TODO) +one_step_pathway ATP_degradation k 0.0000177 /sec Delta_ATP, Singhal et al. Table S2 (TODO) diff --git a/biocrnpyler/utils/plotting.py b/biocrnpyler/utils/plotting.py index 64d5f4ac..2080b861 100644 --- a/biocrnpyler/utils/plotting.py +++ b/biocrnpyler/utils/plotting.py @@ -14,6 +14,8 @@ import statistics from warnings import warn +import numpy as np + from ..components.basic import DNA, RNA, Protein from ..components.dna.cds import CDS from ..components.dna.construct import Construct @@ -21,9 +23,11 @@ from ..components.dna.promoter import Promoter from ..components.dna.rbs import RBS from ..components.dna.terminator import Terminator +from ..core.component import Component from ..core.polymer import OrderedPolymer from ..core.propensities import MassAction from ..core.species import ComplexSpecies, Species +from ..utils.units import mins, mM, nM, uM from . import member_dictionary_search HAVE_MATPLOTLIB = False @@ -1352,3 +1356,261 @@ def render_network_bokeh( raise ValueError("To export you must supply export_name keyword.") export_svgs(plot, filename=export_name + '.svg') return plot + + +def plot_all_species_containing( + results, + crn, + species_list, + time_units=mins, + time_label='min', + concentration_units=uM, + concentration_label='uM', + legend_fontsize='small', + show_material=True, + plot_total=True, + trace_offset=0, +): + """Plot time responses for all species containing a given species. + + Generate a plot showing the concentrations of all species that contain + a given species (or list of species), including complexes. + + Parameters + ---------- + results : DataFrame + Pandas dataframe containing the simulation results. + crn : ChemicalReactionNetwork + CRN being simulutated (used to determine species names). + species_list : List of Species or str + Species to be used to determine full list of species to plot. + time_units : float + Time units factor used to convert from seconds to desired unit. + legend_fontsize : str, default='small' + Font size to use for the legend. + show_material : bool, default=True + Include the material type in the legend string. + plot_total : bool, default=True + Plot the total concentration in addition to individual species. + trace_offset : float, default=0 + Value to use in offseting traces, to avoid overlap. + """ + if not isinstance(species_list, (list, tuple)): + species_list = [species_list] + + timepts = results['time'] / time_units + offset = 0 + for species in species_list: + if isinstance(species, Component): + species = species.species + + compounds = crn.get_all_species_containing(species) + if len(compounds) == 0: + warn(f"No compounds found with species {species}") + return + elif len(compounds) > 1: + species_total = np.zeros(results['time'].size) + for compound in compounds: + plt.plot( + timepts, + results[str(compound)] / concentration_units + offset, + label=compound.pretty_print(show_material=show_material), + ) + species_total += results[str(compound)] + offset += trace_offset + if plot_total: + plt.plot( + timepts, + species_total / concentration_units, + label=f"Total {species.pretty_print()}", + ) + offset += trace_offset + else: + plt.plot( + timepts, + results[str(compounds[0])] / concentration_units, + label=species.pretty_print(), + ) + offset += trace_offset + + plt.legend(fontsize=legend_fontsize) + plt.xlabel(f"Time [{time_label}]") + plt.ylabel(f"Concentration [{concentration_label}]") + + +def plot_gene_expression_data( + results, + crn, + gene_list, + time_units=mins, + time_label='min', + concentration_units=[nM, uM, uM, mM], + concentration_label=['nM', 'uM', 'uM', 'mM'], + legend_fontsize='x-small', + title_fontsize='small', + show_material=True, + plot_total=True, + trace_offset=0, + resource_labels=['metabolite_AAs', 'metabolite_NTPs', 'metabolite_ATP'], +): + """Plot DNA, RNA, and protein concentrations for a list of genes. + + Generate a plot showing the concentrations of all species associated with + expression of a gene, including complex containing the species. + + Parameters + ---------- + results : DataFrame + Pandas dataframe containing the simulation results. + crn : ChemicalReactionNetwork + CRN being simulutated (used to determine species names). + gene_list : List of Species or str + Genes to be used to determine full list of species to plot. + time_units : float, default=bcp.units.mins + Units factor used to plot scale time. + time_label : str, default='min' + String to use in time axis label + concentration_units : float or list of float, default=[nM, uM, uM, mM] + Units factor to scale concentrations of DNA, RNA, protein, resources. + concentration_label : str or list of str, default=['nM', 'uM', 'uM', 'mM'] + Strings to use in labels for DNA, RNA, protein, resources. + legend_fontsize : str, default='x-small' + Font size to use for the legend. + title_fontsize : str, default='small' + Font size to use for individual panel titles. + show_material : bool, default=True + Include the material type in the legend string. + plot_total : bool or list, default=True + Plot the total concentration in addition to individual species. Can + either be True or a list containing 'dna', 'rna', or 'protein'. + trace_offset : float or list of float, default=0 + Value to use in offseting traces, to avoid overlap. If a list is + specified, gives offsets for DNA, RNA, protein, resources. + resource_labels : list of str, optional + Names of resources. Defaults to ['metabolite_AAs', 'metabolite_NTPs', + 'metabolite_ATP']. + """ + if not isinstance(gene_list, (list, tuple)): + gene_list = [gene_list] + + if not isinstance(trace_offset, (list, tuple)): + trace_offset = [trace_offset] * 4 + if not isinstance(concentration_units, (list, tuple)): + concentration_label = [concentration_label] * 4 + if not isinstance(concentration_label, (list, tuple)): + concentration_label = [concentration_label] * 4 + + timepts = results['time'] / time_units + base_offset = 0 + for gene in gene_list: + # DNA + plt.subplot(2, 2, 1) + dna_species = gene.dna + dna_species_list = crn.get_all_species_containing(dna_species) + if len(dna_species_list) > 1: + dna_total = np.zeros(results['time'].size) + offset = base_offset + for species in dna_species_list: + plt.plot( + timepts, + results[str(species)] / concentration_units[0] + offset, + label=species.pretty_print(show_material=False), + ) + dna_total += results[str(species)] + offset += trace_offset[0] + if plot_total is True or 'dna' in plot_total: + plt.plot( + timepts, + dna_total / concentration_units[0], + label=f"Total {dna_species.pretty_print()}", + ) + else: + plt.plot( + timepts, + results[str(gene.dna)] / concentration_units[0], + label=dna_species.pretty_print(), + ) + plt.title("DNA", fontsize=title_fontsize) + plt.ylabel(f"Concentration [{concentration_label[0]}]") + plt.legend(fontsize=legend_fontsize) + + # RNA + plt.subplot(2, 2, 2) + rna_species = gene.transcript + rna_species_list = crn.get_all_species_containing(rna_species) + if len(rna_species_list) > 1: + rna_total = np.zeros(results['time'].size) + offset = base_offset + for species in rna_species_list: + plt.plot( + timepts, + results[str(species)] / concentration_units[1] + offset, + label=species.pretty_print(show_material=False), + ) + rna_total += results[str(species)] + offset += trace_offset[1] + if plot_total is True or 'rna' in plot_total: + plt.plot( + timepts, + rna_total / concentration_units[1], + label=f"Total {rna_species.pretty_print()}", + ) + else: + plt.plot( + timepts, + results[str(gene.transcript)], + label=rna_species.pretty_print(), + ) + plt.title("RNA", fontsize=title_fontsize) + plt.ylabel(f"Concentration [{concentration_label[1]}]") + plt.legend(fontsize=legend_fontsize) + + # Protein + plt.subplot(2, 2, 3) + protein_species = gene.protein + protein_species_list = crn.get_all_species_containing(protein_species) + if len(protein_species_list) > 1: + protein_total = np.zeros(results['time'].size) + offset = base_offset + for species in protein_species_list: + plt.plot( + timepts, + results[str(species)] / concentration_units[2] + offset, + label=species.pretty_print(show_material=False), + ) + protein_total += results[str(species)] + offset += trace_offset[2] + if plot_total is True or 'protein' in plot_total: + plt.plot( + timepts, + protein_total / concentration_units[2], + label=f"Total {protein_species.pretty_print()}", + ) + else: + plt.plot( + timepts, + results[str(protein_species)] / concentration_units[2], + label=protein_species.pretty_print(), + ) + plt.title("Protein", fontsize=title_fontsize) + plt.xlabel("Time [time_label]") + plt.ylabel(f"Concentration [{concentration_label[2]}]") + plt.legend(fontsize=legend_fontsize) + + plt.subplot(2, 2, 4) + offset = 0 + for name, label in zip(resource_labels, ['AAs', 'NTPs', 'ATP']): + plt.plot( + timepts, + results[name] / concentration_units[3] + offset, + label=label, + ) + offset += trace_offset[3] + + plt.title("Resources", fontsize=title_fontsize) + plt.xlabel("Time [time_label]") + plt.ylabel(f"Concentration [{concentration_label[3]}]") + plt.legend(fontsize=legend_fontsize) + + plt.gcf().align_labels() + plt.tight_layout() diff --git a/docs/_autogen_mixtures.rst b/docs/_autogen_mixtures.rst index f66a7cb5..5dcc368e 100644 --- a/docs/_autogen_mixtures.rst +++ b/docs/_autogen_mixtures.rst @@ -49,5 +49,6 @@ Pure :nosignatures: :recursive: + PURE BasicPURE diff --git a/docs/components.rst b/docs/components.rst index 28352682..661f6847 100644 --- a/docs/components.rst +++ b/docs/components.rst @@ -136,4 +136,4 @@ generates its species and reactions during compilation. For example:: pass By defining custom components, users can encode arbitrary logic or kinetic -forms while still integrating seamlessly with mixtures and Mechanisms. +forms while still integrating seamlessly with mixtures and mechanisms. diff --git a/examples/gfp_expression.py b/examples/gfp_expression.py index b96092c4..9c9f0ba3 100644 --- a/examples/gfp_expression.py +++ b/examples/gfp_expression.py @@ -17,8 +17,8 @@ name='gfp', promoter='pconst', rbs='rbs_strong', protein='GFP' ) -# Simulation parmaters -initial_conditions_dict = {'dna_gfp': 1 * nM} +# Simulation parameters +initial_conditions = {'dna_gfp': 1 * nM} timepts = np.linspace(0, 6 * hrs, 1000) # @@ -30,7 +30,7 @@ ) expr_crn = expr_mixture.compile_crn() expr_res = expr_crn.simulate_with_bioscrape_via_sbml( - timepts, initial_condition_dict=initial_conditions_dict + timepts, initial_condition_dict=initial_conditions ) # @@ -42,7 +42,7 @@ ) simple_crn = simple_mixture.compile_crn() simple_res = simple_crn.simulate_with_bioscrape_via_sbml( - timepts, initial_condition_dict=initial_conditions_dict + timepts, initial_condition_dict=initial_conditions ) # @@ -54,7 +54,7 @@ ) regular_crn = regular_mixture.compile_crn() regular_res = regular_crn.simulate_with_bioscrape_via_sbml( - timepts, initial_condition_dict=initial_conditions_dict + timepts, initial_condition_dict=initial_conditions ) # @@ -66,7 +66,7 @@ ) energy_crn = energy_mixture.compile_crn() energy_res = energy_crn.simulate_with_bioscrape_via_sbml( - timepts, initial_condition_dict=initial_conditions_dict + timepts, initial_condition_dict=initial_conditions ) # @@ -144,7 +144,7 @@ plt.figure(4) plt.clf() -cfp_initial_conditions = initial_conditions_dict +cfp_initial_conditions = initial_conditions.copy() cfp_initial_conditions['dna_cfp'] = 1 * nM # Add some additional DNA that will utilize resources @@ -236,7 +236,7 @@ ) pure_crn = pure_mixture.compile_crn() pure_res = pure_crn.simulate_with_bioscrape_via_sbml( - timepts, initial_condition_dict=initial_conditions_dict + timepts, initial_condition_dict=initial_conditions ) plt.figure(5) @@ -248,3 +248,192 @@ plt.xlabel("Time [min]") plt.ylabel("Concentration [uM]") plt.legend() + +# +# PURE mixtures +# + +pure_timepts = np.linspace(0, 180 * min) + +simple_mixture = bcp.PURE( + name='simple', components=[gfp_dna, cfp_dna], + include_machinery=False, include_resources=False, include_fuel=False) +simple_crn = simple_mixture.compile_crn() +simple_res = simple_crn.simulate_with_bioscrape_via_sbml( + pure_timepts, initial_condition_dict=cfp_initial_conditions +) + +machinery_mixture = bcp.PURE( + name='machinery', components=[gfp_dna, cfp_dna], + include_machinery=True, include_resources=False, include_fuel=False) +machinery_crn = machinery_mixture.compile_crn() +machinery_res = machinery_crn.simulate_with_bioscrape_via_sbml( + pure_timepts, initial_condition_dict=cfp_initial_conditions +) + +resource_mixture = bcp.PURE( + name='resources', components=[gfp_dna, cfp_dna], + include_machinery=True, include_resources=True, include_fuel=False) +resource_crn = resource_mixture.compile_crn() +resource_res = resource_crn.simulate_with_bioscrape_via_sbml( + pure_timepts, initial_condition_dict=cfp_initial_conditions +) + +energy_mixture = bcp.PURE( + name='energy', components=[gfp_dna, cfp_dna], + include_machinery=True, include_resources=True, include_fuel=True) +energy_crn = energy_mixture.compile_crn() +energy_res = energy_crn.simulate_with_bioscrape_via_sbml( + pure_timepts, initial_condition_dict=cfp_initial_conditions +) + +energy_mixture_gfp = bcp.PURE( + name='energy', components=[gfp_dna], + include_machinery=True, include_resources=True, include_fuel=True) +energy_crn_gfp = energy_mixture_gfp.compile_crn() +energy_res_gfp = energy_crn_gfp.simulate_with_bioscrape_via_sbml( + pure_timepts, initial_condition_dict=initial_conditions +) + +plt.figure(6) +plt.clf() + +plt.plot( + pure_timepts / min, simple_res['protein_GFP'] / uM, + label='GFP+CFP, simple') +plt.plot( + pure_timepts / min, machinery_res['protein_GFP'] / uM, + label='GFP+CFP, machinery') +plt.plot( + pure_timepts / min, resource_res['protein_GFP'] / uM, + label='GFP+CFP, resources') +plt.plot( + pure_timepts / min, energy_res['protein_GFP'] / uM, + label='GFP+CFP, energy') +plt.plot( + pure_timepts / min, energy_res_gfp['protein_GFP'] / uM, + label='GFP only, energy') +plt.plot( + timepts[0:50] / min, pure_res['protein_GFP'][0:50] + 0.1 / uM, '--', + label='GFP only + 0.1, basic') + +plt.ylim([0, 12]) + +plt.title("PURE Mixture Comparisions") +plt.xlabel("Time [min]") +plt.ylabel("Concentration [uM]") +plt.legend() + +# +# Activators and repressors +# + +TetR = bcp.Protein('TetR') +aTc = bcp.Species('aTc') +TetR_inactive = bcp.ChemicalComplex( + [TetR.species, aTc], mechanisms={'binding': bcp.One_Step_Binding()} +) +ptet_repressed = bcp.RepressiblePromoter('ptet', TetR) +dna_GFP_repressed = bcp.DNAassembly( + 'GFP', promoter=ptet_repressed, rbs='RBS', protein='GFP', length=714 +) + +initial_conditions = {'dna_GFP': 1 * nM, 'protein_TetR': 30 * uM} +repressed_mixture = bcp.PURE( + name='repressed', + components=[dna_GFP_repressed, TetR_inactive], + include_machinery=True, + include_resources=True, + include_fuel=True, +) +repressed_crn = repressed_mixture.compile_crn() + +TetR = bcp.Protein('TetR') +aTc = bcp.Species('aTc') +TetR_inactive = bcp.ChemicalComplex( + [TetR.species, aTc], mechanisms={'binding': bcp.One_Step_Binding()} +) +ptet_regulated = bcp.RegulatedPromoter('ptet', TetR) +dna_GFP_regulated = bcp.DNAassembly( + 'GFP', promoter=ptet_regulated, rbs='RBS', protein='GFP' +) +regulated_parameters = { + ('transcription', 'ptet_leak', 'ktx'): 50, # unbound + ('transcription', 'ptet_TetR', 'ktx'): 0.001, # bound +} + +regulated_mixture = bcp.PURE( + name='regulated', + components=[dna_GFP_regulated, TetR_inactive], + include_machinery=True, + include_resources=True, + include_fuel=True, + parameters=regulated_parameters, +) +regulated_crn = regulated_mixture.compile_crn() + +plt.figure(7) +plt.clf() + +offset = -0.01 +TetR_0 = initial_conditions['protein_TetR'] +for aTc_0 in np.linspace(0, 50*uM, 5): +# aTc_0 = 0 +# for TetR_0 in np.linspace(0, 20*uM, 5): + repressed_res = repressed_crn.simulate_with_bioscrape_via_sbml( + pure_timepts, initial_condition_dict=initial_conditions + | {'aTc': aTc_0} | {'protein_TetR': TetR_0} + ) + lines = plt.plot( + pure_timepts / min, + repressed_res['protein_GFP'] / uM + offset, + label=f"aTc={aTc_0 / uM} uM, TetR={TetR_0 / uM} uM", + ) + + regulated_res = regulated_crn.simulate_with_bioscrape_via_sbml( + pure_timepts, initial_condition_dict=initial_conditions + | {'aTc': aTc_0} | {'protein_TetR': TetR_0} + ) + plt.plot( + pure_timepts / min, + regulated_res['protein_GFP'] / uM + offset, + color=lines[0].get_color(), + linestyle='--', + ) + + offset += 0.05 + +plt.title("Represssed (-) versus Regulated (--)") +plt.xlabel("Time [min]") +plt.ylabel("Concentration [uM]") +plt.legend() + +# +# PURE debugging +# + +repressed_res = repressed_crn.simulate_with_bioscrape_via_sbml( + pure_timepts, initial_condition_dict=initial_conditions + | {'aTc': 37.5 * uM} | {'protein_TetR': 30 * uM} +) + +plt.figure(8) +plt.clf() +bcp.plot_gene_expression_data( + repressed_res, repressed_crn, dna_GFP_repressed, + trace_offset=[0.01, 0.002, 0.1, 0.1]) +plt.suptitle("Gene expression, repressed", fontsize='large') +plt.tight_layout() + +regulated_res = regulated_crn.simulate_with_bioscrape_via_sbml( + pure_timepts, initial_condition_dict=initial_conditions + | {'aTc': 37.5 * uM} | {'protein_TetR': 30 * uM} +) + +plt.figure(9) +plt.clf() +bcp.plot_gene_expression_data( + regulated_res, regulated_crn, dna_GFP_regulated, + trace_offset=[0.01, 0.002, 0.1, 0.1]) +plt.suptitle("Gene expression, regulated", fontsize='large') +plt.tight_layout()