Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Generic stochastic model #1417

Merged
merged 27 commits into from
Jan 19, 2022
Merged

Generic stochastic model #1417

merged 27 commits into from
Jan 19, 2022

Conversation

phumtutum
Copy link
Contributor

@phumtutum phumtutum commented Nov 20, 2021

In this pr:

  • A generic stochastic model is added (MarkovJumpModel) which uses Gillespie's algorithm
  • An example of a model that illustrates how simply other stochastic models can be added from now on (Michaelis Menten model)
  • Added a degradation model for testing as well
  • Removed old degradation model

See also #1067 #890

@codecov
Copy link

codecov bot commented Nov 20, 2021

Codecov Report

Merging #1417 (27453ef) into master (edbf7c1) will not change coverage.
The diff coverage is 100.00%.

Impacted file tree graph

@@            Coverage Diff            @@
##            master     #1417   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           87        90    +3     
  Lines         8899      8926   +27     
=========================================
+ Hits          8899      8926   +27     
Impacted Files Coverage Δ
pints/toy/__init__.py 100.00% <ø> (ø)
pints/toy/stochastic/__init__.py 100.00% <100.00%> (ø)
pints/toy/stochastic/_degradation_model.py 100.00% <100.00%> (ø)
pints/toy/stochastic/_markov_jump_model.py 100.00% <100.00%> (ø)
pints/toy/stochastic/_michaelis_menten_model.py 100.00% <100.00%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update edbf7c1...27453ef. Read the comment docs.

@phumtutum
Copy link
Contributor Author

Now, the old degradation model is removed and the pr is ready for review

@MichaelClerx MichaelClerx mentioned this pull request Nov 30, 2021
Stochastic Toy Problems
***********************

The `stochastic` module provides toy :class:`models<pints.toy.stochastic.MarkovJumpModel>`,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess initially it's just models?

@@ -117,6 +117,7 @@ relevant code.
- [SIR Epidemiology model](./toy/model-sir.ipynb)
- [Stochastic Degradation model](./toy/model-stochastic-degradation.ipynb)
- [Stochastic Logistic model](./toy/model-stochastic-logistic-growth.ipynb)
- [Michaelis Menten model](./toy/model-stochastic-michaelis-menten.ipynb)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Separate category here too?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I think separate out the stochastic ones into their own subsection of "toy models".

@phumtutum
Copy link
Contributor Author

This pr updates the stochastic degradation model. This model is used in a ipynb in #1413 and therefore changes need to be made to that ipynb in this pr once that pr is merged into master. I will make the requested changes at once with those changes.

Copy link
Member

@MichaelClerx MichaelClerx left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fine with me, if @ben18785 likes it!

Copy link
Collaborator

@ben18785 ben18785 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @phumtutum -- looks really good. There are a few things that I think would be good to change, then I'm happy to merge.

In the degradation notebook:

  • I think this line needs rewriting: "averaging the concentration values at each time step, produces a reproducible result" -- I don't think that "reproducible result" is the correct terminology here. I guess you just want to say that the mean of the distribution is that exponential thingy?

  • "The deterministic mean (from above) is plotted with the deterministic standard deviation" -- I think this is a bit ambiguous. I think you want to say, "We now plot the analytic mean and standard deviation of this process."

In the michaelis-menten notebook:

  • "The Michaelis Menten model is a stochastic model, more details can be found here.": the MM model actually isn't necessarily a stochastic model (we have a deterministic version of a similar but not same version of the model in the package), and the stuff you link to is about the deterministic model. So, can we say, "This is a stochastic version of the Michaelis-Menten model, and more details about the deterministic analogue can be found here."
  • Can you also provide a little bit of a description about what this model represents in the first part of this notebook?
  • "This function, given a set of parameters and times, computes the appropriate times and values using Gillispie's algorithm and then uses interpolation to find the values at the exact times requested" -- it's not clear what the values are (i.e. they are counts of molecules) and "Gillespie" has a typo. I think you also need to explain what the interpolation does and why it is needed. I've noticed "Gillespie" is misspelled throughout the notebook actually.
  • "Given the stochastic nature of this model we can use multiple simulations to make sure that the runs are covering the same model with the same parameters. Our simulations are close, suggesting we are obtaining accurate simulations." -- this isn't correct, really. All you are showing is that there is not a large amount of stochasticity in these models for these initial values. Indeed, can we show that when we use much smaller molecule counts as initial conditions that we get much more stochasticity?
  • It would be good to show that the solutions with high molecular counts correspond with the solutions of the corresponding ODE: you could code this up using scipy's odeint within the notebook

@@ -117,6 +117,7 @@ relevant code.
- [SIR Epidemiology model](./toy/model-sir.ipynb)
- [Stochastic Degradation model](./toy/model-stochastic-degradation.ipynb)
- [Stochastic Logistic model](./toy/model-stochastic-logistic-growth.ipynb)
- [Michaelis Menten model](./toy/model-stochastic-michaelis-menten.ipynb)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I think separate out the stochastic ones into their own subsection of "toy models".

if np.any(times < 0):
raise ValueError('Negative times are not allowed.')

variance = np.exp(-2 * k * times) * (-1 + np.exp(k * times)) * self._x0
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I checked -- this is correct

if np.any(times < 0):
raise ValueError('Negative times are not allowed.')

mean = self._x0 * np.exp(-k * times)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Correct

r"""
A general purpose Markov Jump model used for any systems of reactions
that proceed through jumps. We simulate a population of N different species
reacting through M different mechanisms
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"reacting through M different reaction equations"

2. Calculate the time :math:`\tau` until the next single reaction as

.. math::
\tau = \frac{-\ln(r)}{a_0}
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We haven't said what a_0 is.

Comment on lines 34 to 36
3. Decide which reaction, i, takes place using r_1 * a_0 and iterating

through propensities
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This needs quite a bit more explanation, especially since this is the key point of the algorithm

@phumtutum
Copy link
Contributor Author

@ben18785 I have addressed your comments. If there is anything else I should change please let me know.

@ben18785 ben18785 merged commit c7ea751 into master Jan 19, 2022
@MichaelClerx MichaelClerx deleted the i-1067-generic-stochastic-model branch January 19, 2022 15:38
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants