Model Reactions and Implementation
Model overview
The reactions currently incorporated in the model are shown below. The simulator converts all the reactions below into the appropriate rate equations derived from mass action kinetics. This produces a system of ordinary differential equations that is then numerically integrated with the solve_ivp() function in Python’s SciPy package. Concentrations in the simulator are in nM and time is in seconds. The default values for rate constants can be found here
The curly bracket notation specifies the nomenclature of the species where i,j,k,…etc represent the domain numbers with the first index specifying the input domain and the second index specifying the output domain {input,output}. The names of species and the rate constants in the schematics below match how they are defined as variables in the model. The model is matrixed based so the index in the curly brackets also corresponds to the matrix index where the species concentrations are stored, or where the value of individual rate constants are stored in the rate constant matrices. See here for more details on matrix definitions.
The model uses a universal toehold for all strand displacement reactions. So from the model’s perspective all outputs have the same toehold and any output with an output domain index that matches the input domain index of a gate can react with that gate, i.e., O{3,1} will react with G{1,4} and G{1,5}. If it is desired to model a system in which multiple toeholds with different strand displacement rates, the appropriate index of the krsd rate constant matrix can be changed in molecular_species for each gate.
The model treats input species as outputs with identical input-output indices, i.e., I{1} = O{1,1}, I{4} = O{4,4}. The user can specify input species in the model functions, such as when defining species in molecular_species() or when pulling concentrations for analysis with output_concentration(), but internally these will be handled in the output matrices. This implementation means other species that produce outputs with identical input-output indices, such as G{1,1}, AG{3.2,2} should not be used in the model as their outputs will behave as inputs rather than outputs.
There are a few things to note about the way the rate constant matrices are defined and implemented.
For RNA strand displacement reactions, unless otherwise stated, the forward rate constants are associated with the gates. This means the rate constant for a forward strand displacement reaction for a specific gate, for example G{1,3}, can be changed and that rate constant will be used with all outputs that react with that gate, i.e., O{4,1} and O{5,1} will both use the forward rate constant associated with G{1,3}. It is not possible to give O{4,1} and O{5,1} unique forward rate constants for a strand displacement reaction with G{1,3} (or any gate with input domain index 1). Changes to forward RNA strand displacement reactions for individual gates should be done in the molecular_species() function when defining the gates. Changing the forward RNA strand displacement rate constant is not a valid option for outputs in molecular_species().
The reverse rate constants for RNA strand displacement reactions are associated with the outputs. This means the rate constant for a reverse strand displacement reaction for a specific output, for example O{1,3}, can be changed and that rate constant will be used with all GO{} species that react with the output, i.e., GO{4,1} and GO{5,1} will both use the reverse rate constant associated with O{1,3}. It is not possible to give GO{4,1} and GO{5,1} unique reverse rate constants for a strand displacement reaction with O{1,3} (or any output with input domain index 1). Changes to reverse RNA strand displacement reactions for individual outputs should be done in the molecular_species() function when defining the outputs, or when defining the gates – the reverse rate that is defined with a gate with correspond to the output of that gate. NOTE: because inputs are stored on the diagonal of the output matrix the diagonal of the reverse RNA strand displacement rate constant matrix is set to 0. This is because inputs only possess an output domain index and cannot participate in reverse RNA strand displacement reactions
ctRSD Reaction Schematics
Below the subscript “temp” indicates a DNA transcription template that encodes for the RNA species specified in its name.
We model RNA transcription and degradation using a first order approximation for enzyme kinetics. From this approximation the transcription rate is linearly proportional to the template concentration, e.g., ktxnO{k,i}*[Otemp{k,i}] where ktxn{k,i} is the apparent first order rate constant for production of RNA from Otemp{k,i}. Likewise, the degradation rates (if included in the model) are linearly proportional to the concentration of RNA, e.g., kssdO{i,j}*[O{i,j}].