Help with Model

Question

[Stele]

Hey everyone, I am coding a stochastic model in python and I have hit a wall and don't know what to do or how to code something.

This model is to stochastically simulate the interactions between 2 molecules, molecule1 and molecule 2, taking into account their transcription, splicing, and translation rates (and degradation). It is meant to have an ON/OFF switch that activates/deactivates the transcription of the genes in question. My issue is with the probability of the switch being ON. I want it to have a base rate and increase with the amount of protein present (say, increases 10% every 10 proteins with a cap on 30%). There is a loop of positive autofeedback for each molecule, as well as cross-feedback.

If anyone has any idea on how to code for this I would appreciate a comment. I can share the full code if requested (don't want to share the complete thing here, as it would make the post massive).

Thanks a lot!

[Joe]

I would imagine this is something one of the AI code tools could handle.

Answer 1

[i.sudbery]

I'm not quite sure what you are asking. Stochastic reaction networks like this can be coded using the Gillespie algorithm.

Genes being "on" or "off" is usually imaged as two different states, with the transition between the two dependent on the regulator.

In the simplest case you would have a model where a regulator binds to a gene, which switches the gene on, which can then transcribe a primary RNA, which can be degraded, or can be spliced to make mRNA, which can be degraded or translated to make Protein (which doesn't consume the mRNA). The Protein can also be degraded.

Reaction                        Propensity
Reg + Gene_off -> Gene_on       k_on * [Reg] * [Gen_off]
Gene_on -> Reg + Gene_off       k_off * [Gene_off]
Gene_on -> Gene_on + RNA        k_transcription * [Gene_on]
RNA -> 0                        delta_rna * [RNA]
RNA -> mRNA                     k_splicing * [RNA]
mRNA -> 0                       delta_mRNA * [mRNA]
mRNA -> mRNA + Protein          k_translation * [mRNA]
Protein -> 0                    delta_protein * [Protein]

You'd have to figure out the parameter values for k_on, k_off and [Reg] that gave you the properties you wanted. And of course, this system response linearly to regulator protein, you'd need a more complex system if you wanted to the reactions to be non-linear.

To simulate under Gillespie you use total propensity to calculate a waiting time for the next event based on an poisson waiting time model, and then a draw from a multinomial uniform distribution, weighted by the propensities to decided what reaction happens, update entity counts, and repeat.

I have to say, this does sound awefully like something someone might set as homework for a grad level systems biology course. In fact, I seem to remember having to model something very similar for SB100 taught by Johan Paulson.