This example shows how the Repressilator model can be used.
The model, formulated as an ODE, has 6 state variables: 3 mRNA concentrations and 3 protein concentrations. Only the mRNA concentrations are visible. In the example below we'll call the three outputs m-lacI, m-tetR, and m-cl.
The model has 4 parameters, alpha_0, alpha, beta, and n.
For an analysis using ABC, see Toni et al..
In [1]:
import pints
import pints.toy
import pints.plot
import numpy as np
import matplotlib.pyplot as plt
# Create a model
model = pints.toy.RepressilatorModel()
# Run a simulation
parameters = model.suggested_parameters()
times = model.suggested_times()
values = model.simulate(parameters, times)
print('Parameters:')
print(parameters)
# Plot the results
plt.figure()
plt.xlabel('Time')
plt.ylabel('Concentration')
plt.plot(times, values)
plt.legend(['m-lacI', 'm-tetR', 'm-cl'])
plt.show()
With these parameters, the model creates wide AP waveforms that are more reminiscent of muscle cells than neurons.
We now set up a simple optimisation problem with the model.
In [2]:
# First add some noise
sigma = 5
noisy = values + np.random.normal(0, sigma, values.shape)
# Plot the results
plt.figure()
plt.xlabel('Time')
plt.ylabel('Concentration')
plt.plot(times, noisy)
plt.show()
Next, we set up a problem. Because this model has three outputs, we use a MultiOutputProblem.
In [3]:
problem = pints.MultiOutputProblem(model, times, noisy)
loglikelihood = pints.GaussianKnownSigmaLogLikelihood(problem, sigma)
Now we're ready to try some inference, for example MCMC:
In [4]:
# Initial guesses
x0 = [
[2, 800, 3, 3],
[1, 1200, 6, 1],
[3, 2000, 1, 4],
]
mcmc = pints.MCMCController(loglikelihood, 3, x0)
mcmc.set_log_to_screen(False)
mcmc.set_max_iterations(6000)
chains = mcmc.run()
We can use rhat criterion and look at the trace plot to see if it looks like the chains have converged:
In [5]:
# Check convergence using rhat criterion
print('R-hat:')
print(pints.rhat_all_params(chains))
plt.figure()
pints.plot.trace(chains, ref_parameters=parameters)
plt.show()
So it seems MCMC gets there in the end!
We can use the final 1000 samples to look at the predicted plots:
In [6]:
samples = chains[1][-1000:]
plt.figure(figsize=(12, 6))
pints.plot.series(samples, problem)
plt.show()