Notebook to infer reaction constants of a stochastic enzyme reaction

[SidRama]

Summary:

Notebook demonstrating the use of ELFI to perform ABC inference of reactions constants of a toy stochastic enzyme reaction.

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company):

By submitting this pull request, the copyright holder is agreeing to license the submitted work under the following licenses:

• Code: BSD3 (https://opensource.org/licenses/BSD-3-Clause)
• Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

[OssiGalkin]

Random_state argument is used in definition of simulator, but it is no actually used, therefore result will vary between runs.

running notebook is incredible slow because batch_size=1. I didn't even had patience to wait it complete ;) I don't blame anyone but Elfi for this, as using batch_size over one is too hard at the moment.
Non-Aalto students probably don’t know what is Triton, also it would nice to see code that does Triton part. For example, code for this part would be interesting to see:
“Similar to rejection sampling, we can also play around with the Summary statistics. In the next analysis, the summary statistics with lag greater than 1 were removed. That is, we omit S2, S4, S6 and S8. We draw 4000 samples per chain and again run it on Triton as it is time consuming. Below, we simply load the samples.”

Different autocovariance functions could be defined in single function:

# Summary statistic for P
def autocov(x, lag=1, i=0):
    C = np.mean(x[i,lag:] * x[i,:-lag])
    return C

and then give correct value for i, for example:
S5 = elfi.Summary(autocov_3, X, name='S5')
would become:
S5 = elfi.Summary(autocov, X, i=3, name='S5')

[cagatayyildiz]

I think this notebook is prepared really well. Below are my comments:

• I wonder whether the priors oversimplify the problem or not. What are the realistic ranges for c1, c2 and c3 (if there is any)?

• The priors specify the range for each varible as [0,2]. Yet, the inferred values go from 0 to 3. Why is so?

• Is there a reason why the time in the plots? The code and the plot in "True model of the stochastic enzyme reaction" plot seemingly do not match.

• Isn't threshold printed by [20] suspiciously high? (see the second cell in Inference with rejection sampling section)

• The number of samples in [25] is still 1000, not 100000 as written in the cell above. Is it the number of simulations (rather than the number of samples) that is wanted to be set to 1000?

[vuolleko]

Excellent and very well documented notebook!

Some comments for consideration:

• In gillespie, there's unnecessary mix of +-1 with index
• Isn't the convention usually that rates (and their priors) are in [0,1] and the stoichiometry matrices handle splitting/merging of molecules?
• c3 could be a hierarchical prior that depends on c2 (since c2+c3 <= 1)
• For some reason the labelling in figure resulting from In [6] does not match the code???
• Same figure^, maybe comment why the reactions happen very fast during the first ms
• You could use 1 autocov function for all 4 species
• "Note: number of samples = N/quantile." Perhaps better to use the word "proposals" to avoid confusion between samples from prior vs. posterior
• Priors were set in range [0,2], so why plots from In [21] have range [0,3]?
• What was the choice of summary statistics based on? The output from simulator is not particularly noisy, so calculating distance direct from output might work as well.
• Awesome, I appreciate the use of the Triton cluster, but the notebook should preferably be self-contained; it's fine to mention results obtained elsewhere, but loading non-existent files is maybe not
• Move to zoo/notebooks/enzyme_reaction.ipynb or something

[vuolleko]

@OssiGalkin It's hardly the fault of ELFI if the code is slow with batch_size=1, but then again if you are able to use a larger batch_size, the speedup offered by ELFI should be appreciated.

See code in elfi.examples.lotka_volterra for a somewhat similar case in vectorized form.

[SidRama]

Thank you everyone for the comments! I will make the updates soon.