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.

[Stele]

Hey, thanks for the reply! I think I didn't make the question very understandable. I do know how to code for this, the wall I have hit is "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%)."

[i.sudbery]

There is no way to "code" this directly. Your problem is mathematics, not programming.

Standard stochastic models deal with a set of states, and propensities to transition between them. The propensities are generally modelled as the product of count of each molecule involved and a rate constant (assuming mass action kinetics). If you want your system to exhibit a particular property, you need to fit parameters to the model until it exhibits that property, or alter the set of possible states and transitions until it does what you want it to do.

Its unclear what you mean in your question by a "switch". Do you mean that a gene is on? Or are you talking about the system occupying one of two possible stable states, one of which is labelled "on" (e.g. a bistable negative feedback loop). I guess what I'm asking,is what is the physical nature of your "switch"?

Assuming the first (although the following also applies to the second), in any given cell, a gene is either on or off - probabilities apply to populations, not to individuals. So we can say that the gene is on in 30% of cells in any given instant, or we can say the gene is on 30% of the time, but we can't say that there is a 30% chance of a gene being on in a specific cell in a given instant. Even for a gene, being "ON" could mean two different things. We could mean that in 30% of model runs the gene is in the ON state at a given time point, or it could mean that when the gene is in the ON state, there is a 30% chance that it generates a transcript in some time window.

I can think of no mechanistic system that would give you the behavoir you want if we are talking about the gene being in the ON state.

The simpliest system for a gene moving between an on and an off state is:

Reaction                        Propensity
Reg + Gene_off -> Gene_on       k_on * [Reg] * [Gen_off]
Gene_on -> Reg + Gene_off       k_off * [Gene_on]

The long term average of the number of gene copies in the on state here is k_on*[Reg]/k_off (I think, although I'd check that). But there is no way to cap that at 30% of gene copies being "ON" or 30% of cells have genes in the ON state.

If by 30% you mean that there is a 30% probability of a transcript being produced in a second (for example), you could alter k_transcription until it gave the required probability (i guess if you set k_transcription to 0.3/s then you'd get a 30% of firing within a 1 second period).