8.7. Developer Overview

Introduction

BioCRNpyler is a software tool designed to rapidly compile large biological chemical reaction networks (CRNs) from simple user specifications (written in python). It has built in support of a number of models for transcription, translation, and gene expression regulation using components common in E. coli synthetic biology. This tutorial explains the inner workings of BioCRNpyler and shows how to create custom mixtures, components, and mechanisms. Specifically, we will go through making a custom gene expression model: >\(G \to G + P \rightleftharpoons G:P \to G + P + X\)

here \(G\) is a gene and \(P\) is a polymerase and \(X\) is the protein expressed by \(G\). No translational machinery is included in this model, making it one of the simplest possible for expression. Note that we are ignoring translation for simplicitiy, not becuase it isn’t important.

On the top level, BioCRNpyler uses three kinds of objects:

  • Mechanisms: are the details of how a physics process is implemented as a CRN. These take the form of black box reaction schemas which compile into a CRN containing all the intermediate steps required to get from a specified input to an output.

  • Components: are the ingredients one might imagine adding to a test tube, say from a pipette. They do not include all chemical species involved in a reaction, but just the key ones we might experimentally modulate. Components may contain their own mechanisms or default to those used by a mixture. An example of a component is a piece of DNA encoding a gene. A DNA-Transcription factor complex, on the other hand, would not normally be a component.

  • Mixtures: can be thought of as the “reaction soup” we are working in. Mixtures contain default components and mechanisms. Components are added to mixtures to create different reaction conditions.

Internally, BioCRNpyler tells the Mixture to compile all its Components. Each Component contains its own Mechanisms (or defaults to Mechanisms defined in the Mixture) and calls each Mechanism (read: reaction schema) to generate a set of chemical species and reactions which are combined into a complete CRN. BioCRNpyler also has its own internal CRN representation, which we will discuss next.

Chemical Reaction Network (CRN) model

A CRN is a set of species \(S\) and a set of reactions \(R\) where each reaction is expressed \(I \rightarrow O\) where \(I\) are the inputs species, \(O\) are the output species. Each reaction occurs with a rate function (propensity) \(\rho(x)\) which takes the state of the CRN (the values of all the species) as an input. By default, reactions use massaction rates: \(\rho(x) = k \Pi_{s \in I} x_s\) here \(k\) is some constant and \(x_s\) is the value of the species \(s\). A number of built in propensities exist and are described in the documentation, including a general propensity allowing for an arbitrary function.

Internally, BioCRNpyler represents species as strings involving a type identifier and a name: type_name. This is to allow for species to be identified as “gene_X”, “mrna_X”, etc. Complexes between species can be created automatically using the ComplexSpecies constructor or given custom defined names. By default, a complex of gene_X and protein_Y would be called complex_gene_X_protein_Y. This would be considered different from complex_protein_Y_gene_X in Bioscrape’s CRN semantics because species here are effectively strings.

Reactions are stored as lists of species (for the inputs and outputs) and a rate constant k. Non massaction reactions also require a parameter dictionary of their relevant parameter values. Massaction reactions are allowed to be reversible, in which case they are thought of as two irreversible reactions. Reaction rates default to 1.0.

Now, we will create the CRN described above directly and approximate it with a non-massaction propensity.

[1]:

from biocrnpyler.core.chemical_reaction_network import ChemicalReactionNetwork
from biocrnpyler.core.species import Species, Complex
from biocrnpyler.core.reaction import Reaction
from biocrnpyler.core.propensities import HillPositive
#create the three species in the CRN
G = Species(name = "G", material_type = "dna")
P = Species(name = "P", material_type = "protein")
X = Species(name = "X", material_type = "protein")
PG = Complex([P, G]) #complex takes a list of species and returns a complex
species = [P, G, X, PG] #a list of species

#Create the reversible reaction: + P <--> G:P
kf = 100 #Forward reaction rate
kr = .01
inputs1 = [G, P]
outputs1 = [PG]
rxn1 = Reaction.from_massaction(inputs1, outputs1, k_forward = kf, k_reverse = kr) #type defaults to massaction
#Create the irreversible reaction G:P --> G + P + X
inputs2 = [PG]
outputs2 = [G, P, X]
kexpress = 1.
rxn2 = Reaction.from_massaction(inputs2, outputs2, k_forward = kexpress)



rxns = [rxn1, rxn2] #a list of reactions

CRN = ChemicalReactionNetwork(species, rxns)

#Species, reactions, and CRNs can all be directly printed
print("species representation:\n", [a.material_type for a in species])
print("\nrxns representation:\n", rxns)
print("\nCRN Representation:\n", CRN)

#We will now create a third reaction which models the production of X as a positive hill function of P
inputs3 = [G, P]
outputs3 = []
khill = 10
#parmeters can be numbers or strings
pos_hill = HillPositive(k=khill, K=0.3, n=2, s1=P)
rxn3 = Reaction(inputs3, outputs3,propensity_type=pos_hill)
CRN2 = ChemicalReactionNetwork(species, [rxn3])
print("\nCRN2:\n",CRN2)

species representation:
 ['protein', 'dna', 'protein', 'complex']

rxns representation:
 [dna[G]+protein[P] <--> complex[dna[G]:protein[P]], complex[dna[G]:protein[P]] --> dna[G]+protein[P]+protein[X]]

CRN Representation:
 Species = protein_P, dna_G, protein_X, complex_dna_G_protein_P_
Reactions = [
        dna[G]+protein[P] <--> complex[dna[G]:protein[P]]
        complex[dna[G]:protein[P]] --> dna[G]+protein[P]+protein[X]
]

CRN2:
 Species = protein_P, dna_G, protein_X, complex_dna_G_protein_P_
Reactions = [
        dna[G]+protein[P] -->
]

Creating a Custom Mechanism: GeneExpression

To create custom Mechanism objects, subclass the Mechanism class and rewrite the object constructor, the update_species function, and the update_reactions function. Briefly:

  • In the constructor we will set the name of the mechanism and the name of the polymerase species, rnap.

  • In update_species, we will create a list of all the species used in the reaction schema: the gene, gene-rnap complex, and the product species.

  • In update_reactions we create a list of all the reactions required for our reaction schema: the polymerase binding and unbinding reactions as well as the reaction producing the gene product X.

Note that this code could be generated much faster using the built in MichaelisMentenRXN Mechanism, but we will do it by hand here for educational purposes.

[2]:

from biocrnpyler.core.mechanism import Mechanism



class GeneExpression(Mechanism):
    #Overwrite the constructor.
    #    Name: the name of the Mechanism (set when it is instantiated).
    #    rnap: the polymerase, which we will allow to be multiple types of object for user convenience
    #    type: this is the "kind" of mechanism - used as a key in mechanism dictionaries
    def __init__(self, name, rnap, type = "gene_expression", **keywords):
        #Check if the rnap type species (see chemical reaction network details below)
        if isinstance(rnap, Species):
            self.rnap = rnap
        else:
            raise ValueError("'rnap' must be a Species!")
        #The superclass constructor will take care of the name
        Mechanism.__init__(self = self, name = name, mechanism_type = type, **keywords) #MUST CALL THE SUPER CONSTRUCTOR!

    #Overwrite update_species:
    #    dna: the name of the gene to be expressed
    #    product: the name of the gene product
    #update_species returns a list of all species used by the reaction schema
    def update_species(self, dna, product):
        #We do not need to do a check on the DNA or product types because that will be performed at the Component level.
        #Create the list of species to return
        species = [dna, self.rnap, product]
        #The Complex returns a ComplexSpecies made up a list of species
        species += [Complex([dna, self.rnap])]
        #Return a list of species
        return species

    #Overwrite update_species:
    #    dna: the name of the gene to be expressed
    #    product: the name of the gene product
    #    component and part_id are used for the mechanism to find parameters approrpiately
    #update_species returns a list of all species used by the reaction schema
    #update_reactions will require rates as well as the relevant species. Returns a list of chemical_reaction_network.reaction
    def update_reactions(self, dna, product, component, part_id = None):

        #Component.get_parameter will automatically search the ParameterDatabases for the best parameter to use.
        #The string names here, 'kexpress', 'kb', 'ku', must be defined by you to match the parameter data file.
        #see parameter jupyter notebook for more information.
        kexpress = component.get_parameter("kexpress", part_id = part_id, mechanism = self)
        kb = component.get_parameter("kb", part_id = part_id, mechanism = self)
        ku = component.get_parameter("ku", part_id = part_id, mechanism = self)

        #complex specie
        comp = Complex([dna, self.rnap])
        #Binding Reaction: dna + rnap <--> dna:rnap
        binding_rxn = Reaction.from_massaction(inputs=[dna, self.rnap], outputs=[comp], k_forward=kb, k_reverse=ku)
        #Catalytic Reaction: dna:rnap --> dna + rnap + product
        cat_rxn = Reaction.from_massaction(inputs=[comp], outputs=[dna, product, self.rnap], k_forward=kexpress)
        #Return a list of reactions
        return [binding_rxn, cat_rxn]

Creating a Custom Component: Gene

To create custom Component objects, subclass the Component class and rewrite constructor, update_species, and update_reactions functions.

  • The Constructor: will set the name of the DNA specie and the name of the protein product

  • update_species: will call each mechanism (in this case just GeneExpression) to get their species

  • update_reactions: will call each mechanism (in this case just GeneExpression) to get their reactions

In general, each component’s functions update_species and update_reactions need to know (via you, the programmer) what the names of the mechanisms they are expected to use are. These mechanisms will be automatically inherited from the Mixture object the Component is added to (by default) but can also be overwritten with the mechanisms keyword in the Component constructor.

[3]:

from biocrnpyler.core.component import Component

class Gene(Component):
    #OVERWRITE CONSTRUCTOR
    def __init__(self, dna, product = None, **keywords):
        #check types for name and product and set internal variables
        #self.internal_species = Component.set_species(species, material_type = None, attributes = None)
        #is a helper function that allows for species to be strings, Species, or Components.
        self.dna = self.set_species(dna, material_type = "dna")

        if product is None: #provide default name for the product
            self.product = self.set_species(self.dna.name, material_type = "protein")
        else:
            self.product = self.set_species(product)

        Component.__init__(self = self, name = dna, **keywords) #MUST CALL THE SUPERCLASS CONSTRUCTOR!


    #OVERWRITE update_species
    def update_species(self):
        #The Component will automatically search for a mechanism called "gene_expression", which it can find in 2 ways
        #    1: it can inherit this from its Mixture (which requires the Mixture has an appropriate "gene_expression" mechanism)
        #    2: this can be passed into the Gene constructor in a dictionary as a keyword arg mechanisms= {'gene_expression':Mechanism [Object Instance]}
        mech_express = self.get_mechanism("gene_expression") #argument is the mechanism type

        #Return the species from the mechanisms in your mixture. In this case, just one.
        return mech_express.update_species(self.dna, self.product)

    #OVERWRITE update_reactions
    def update_reactions(self):

        #argument is the mechanism type
        mech_express = self.get_mechanism("gene_expression")

        #Return the reactions from each mechanism in your mixture. In this case, just this one.
        return mech_express.update_reactions(self.dna, self.product, component = self, part_id = self.name)

Creating a Custom Mixture: ExpressionMixture

To create custom Mixture objects, subclass the Mixture class and rewrite the object constructor function to contain the appropriate default mechanisms and components. All other functionalities will be inherited from the Mixture super class.

[4]:

#ExpressionMixture
from biocrnpyler import Mixture, Protein

class ExpressionMixture(Mixture):
    #OVERWRITE THIS METHOD
    def __init__(self, name="", rnap = "RNAP", **keywords):
        #MUST CALL THE SUPERCLASS CONSTRUCTOR!
        Mixture.__init__(self, name = name, **keywords)

        #RNAP is a component which will be added to the Mixture
        self.rnap = Protein(rnap)

        #add the components to the Mixture
        self.add_components([self.rnap])

        #Create an instance of the GeneExpression mechanism
        mech_express = GeneExpression("gene_expression", self.rnap.get_species())
        #notice the Species inside the Component is passed in with Component.get_species()

        #Create default mechanism dict
        default_mechanisms = {
            mech_express.mechanism_type:mech_express
        }

        #add the mechanisms to the Mixture
        self.add_mechanisms(default_mechanisms)

Combine everything and compile a CRN and print it.

[5]:

#Create a fake parameter dictionary for the example
parameters = {("gene_expression","Reporter", "kexpress"):1.0,
              ("gene_expression","Reporter", "ku"):.01,
              ("gene_expression","Reporter", "kb"):100.0 }

#Instantiate a gene
G1 = Gene("Reporter", "GFP", parameters = parameters)
myMixture = ExpressionMixture(components = [G1])

print(myMixture)

CRN = myMixture.compile_crn()
#Print the CRN
print("Internal String Representation of the CRN:\n", CRN)

print("\nFancier Pretty Print Representation Can also be used for Species, Reactions, and CRNS:\n",
      CRN.pretty_print(show_rates = True, show_material = True, show_attributes = True))

ExpressionMixture:
Internal String Representation of the CRN:
 Species = dna_Reporter, protein_RNAP, GFP, complex_dna_Reporter_protein_RNAP_
Reactions = [
        dna[Reporter]+protein[RNAP] <--> complex[dna[Reporter]:protein[RNAP]]
        complex[dna[Reporter]:protein[RNAP]] --> dna[Reporter]+GFP+protein[RNAP]
]

Fancier Pretty Print Representation Can also be used for Species, Reactions, and CRNS:
 Species(N = 4) = {
    complex[dna[Reporter]:protein[RNAP]] (@ 0),
    dna[Reporter] (@ 0),
    protein[RNAP] (@ 0),
    GFP (@ 0),
}

Reactions (2) = [
0. dna[Reporter]+protein[RNAP] <--> complex[dna[Reporter]:protein[RNAP]]
 Kf=k_forward * dna_Reporter * protein_RNAP
 Kr=k_reverse * complex_dna_Reporter_protein_RNAP_
  k_forward=100.0
  found_key=(mech=gene_expression, partid=Reporter, name=kb).
  search_key=(mech=gene_expression, partid=Reporter, name=kb).
  k_reverse=0.01
  found_key=(mech=gene_expression, partid=Reporter, name=ku).
  search_key=(mech=gene_expression, partid=Reporter, name=ku).

1. complex[dna[Reporter]:protein[RNAP]] --> dna[Reporter]+GFP+protein[RNAP]
 Kf=k_forward * complex_dna_Reporter_protein_RNAP_
  k_forward=1.0
  found_key=(mech=gene_expression, partid=Reporter, name=kexpress).
  search_key=(mech=gene_expression, partid=Reporter, name=kexpress).

]

[6]:

# End