<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 HES1 Michaelis-Menten toy model
The Hes1Model describes the expression level of the transcription factor Hes1, details of the model see here.
The model has 3 state variables $m$, $p_1$ and $p_2$, where $m$ represents the molecular concentrations of Hes1 mRNA, and $p_1, p_2$ represent the molecular concentrations of cytoplasmic and nuclear proteins, respectively.
It is described by the ODE
Step1:  We can get some suggested parameters and plot some data! In experiments, typically it is not easy to have high sampling rate (time points), so our suggested time points match what is feasible to see in actual experiments. However, it is also good to have a look at what the model does by giving it a higher sampling rate (time points).
Step2:  In this set-up, we only see one of the state variables representing the Hes1 concentration. However, it is also interesting to see the dynamics of the whole system (i.e. by inspecting all state variables). This can be done by simulate_all_states(parameters, times) provided in this model.
Step3:  We can explore the model further by varying the initial conditions
Step4:  Then by varying the parameters
Step5:  We can then try to set up and run a MCMC routine with suggested values which origianlly obtained from real experimental results (taken from [1])
Step6:  We can also compare the predictions with these values to what we found
Step7:  This is interesting, as this shows that with the sparse sampling points, we cannot retreive the same dynamical behaviour as the literature model. Our model predictions from the MCMC posterior actually give a closer alignment to the data points.
Since this model has sensitivities, we can try Hamiltonian Monte Carlo for inference.
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
import pints
import pints.toy
model = pints.toy.Hes1Model()
print('Outputs: ' + str(model.n_outputs()))
print('Parameters: ' + str(model.n_parameters()))
Explanation: HES1 Michaelis-Menten toy model
The Hes1Model describes the expression level of the transcription factor Hes1, details of the model see here.
The model has 3 state variables $m$, $p_1$ and $p_2$, where $m$ represents the molecular concentrations of Hes1 mRNA, and $p_1, p_2$ represent the molecular concentrations of cytoplasmic and nuclear proteins, respectively.
It is described by the ODE:
\begin{align}
        \frac{dm}{dt} &= -k_{deg}m + \frac{1}{1 + (p_2/P_0)^h} \
        \frac{dp_1}{dt} &= -k_{deg} p_1 + \nu m - k_1 p_1 \
        \frac{dp_2}{dt} &= -k_{deg} p_2 + k_1 p_1
\end{align}
where $k_1$ is the rate of transport of Hes1 protein into the nucleus; $P_0$ is the amount of Hes1 protein in the nucleus, when the rate of transcription of Hes1 mRNA is at half its maximal value; $\nu$ is the rate of translation of Hes1 mRNA, and $h$ is the Hill coefficient. They are the model parameters to be inferred; whereas $k_{deg}$ is the Hes1 protein degradation rate which is assumed to be known. 
In this model, only the state variable $m$ is an observable, i.e. can be observed in the experiments. The input order of parameters of interest is ${ P_0, \nu, k_1, h }$.
End of explanation
times = model.suggested_times()
smooth_times = np.linspace(times[0], times[-1], 1000)
parameters = model.suggested_parameters()
print(parameters)
plt.figure()
plt.xlabel('Time [minute]')
plt.ylabel('Hes1 concentration')
plt.plot(times, model.simulate(parameters, times), 'o', label='Sparse sampling')
plt.plot(smooth_times, model.simulate(parameters, smooth_times), '--', label='Underlying model')
plt.legend()
plt.show()
Explanation: We can get some suggested parameters and plot some data! In experiments, typically it is not easy to have high sampling rate (time points), so our suggested time points match what is feasible to see in actual experiments. However, it is also good to have a look at what the model does by giving it a higher sampling rate (time points).
End of explanation
all_states = model.simulate_all_states(parameters, smooth_times)
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel('m')
ax.set_ylabel('p1')
ax.set_zlabel('p2')
plt.plot(all_states[:, 0], all_states[:, 1], all_states[:, 2])
plt.show()
Explanation: In this set-up, we only see one of the state variables representing the Hes1 concentration. However, it is also interesting to see the dynamics of the whole system (i.e. by inspecting all state variables). This can be done by simulate_all_states(parameters, times) provided in this model.
End of explanation
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel('m')
ax.set_ylabel('p1')
ax.set_zlabel('p2')
suggested_init = np.array(model.initial_conditions())
suggested_hidden_param = np.array(model.fixed_parameters())
for x in [1.5, 1.6, 2, 2.5, 3, 4, 5]:
    model.set_initial_conditions(suggested_init * x)
    model.set_fixed_parameters(list(suggested_hidden_param[0:-1] / x) + [suggested_hidden_param[-1]])
    all_states = model.simulate_all_states(parameters, smooth_times)
    plt.plot(all_states[:, 0], all_states[:, 1], all_states[:, 2])
plt.show()
# Reset initial conditions
model.set_initial_conditions(suggested_init)
model.set_fixed_parameters(suggested_hidden_param)
Explanation: We can explore the model further by varying the initial conditions:
End of explanation
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel('m')
ax.set_ylabel('p1')
ax.set_zlabel('p2')
# Try some different parameter sets {P0, v, k1, h}
Try_params = [[0.5, 0.5, 0.6, 4.0],
              [0.5, 0.1, 0.01, 5.5],
              [0.4, 0.01, 0.01, 6.0],
              [0.05, 0.05, 0.01, 7.0]]
# Compare with suggested parameters
all_states = model.simulate_all_states(parameters, smooth_times)
plt.plot(all_states[:, 0], all_states[:, 1], all_states[:, 2], label='Suggested parameters')
for i, try_param in enumerate(Try_params):
    all_states = model.simulate_all_states(try_param, smooth_times)
    plt.plot(all_states[:, 0], all_states[:, 1], all_states[:, 2], label='Trial parameters %d'%(i+1))
plt.legend()
plt.show()
Explanation: Then by varying the parameters:
End of explanation
# Get suggested values
values = model.suggested_values()
plt.figure()
plt.xlabel('Time [minute]')
plt.ylabel('Hes1 concentration')
plt.plot(smooth_times, model.simulate(parameters, smooth_times), label='Literature model')
plt.plot(times, values, '*', label='Suggested values')
plt.legend()
plt.show()
# Create an object with links to the model and time series
problem = pints.SingleOutputProblem(model, times, values)
# Create a log posterior
lower = list(parameters * 0.1) + [0.1]
upper = list(parameters * 10) + [5]
log_prior = pints.UniformLogPrior(lower, upper)
log_likelihood = pints.GaussianLogLikelihood(problem)
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
# Run MCMC on the noisy data
x0 = [list(parameters*1.2) + [0.5]]*3
mcmc = pints.MCMCController(log_posterior, 3, x0)
mcmc.set_max_iterations(10000)
chains = mcmc.run()
results = pints.MCMCSummary(
    chains=chains,
    time=mcmc.time(),
    parameter_names=["P0", "nu", "k_1", "h", "sigma"])
print(results)
import pints.plot
pints.plot.trace(chains, ref_parameters=list(parameters)+[0.75])
plt.show()
Explanation: We can then try to set up and run a MCMC routine with suggested values which origianlly obtained from real experimental results (taken from [1]):
[1] Silk, D. et al. Designing attractive models via automated identification of chaotic and oscillatory dynamical regimes. Nat. Commun. 2:489 doi: 10.1038/ncomms1496 (2011).
End of explanation
# Select first chain
chain1 = chains[0]
# Remove burn-in
chain1 = chain1[5000:]
# Plot some predictions with these samples
plt.figure()
plt.xlabel('Time [minute]')
plt.ylabel('Hes1 concentration')
for posterior_param in chain1[::25]:
    model_prediction = model.simulate(posterior_param[:-1], smooth_times)
    plt.plot(smooth_times, model_prediction, c='Gray', alpha=0.3)
plt.plot(smooth_times, model_prediction, c='Gray', alpha=0.3, label='Model prediction')
plt.plot(times, values, 'kx', label='Original data')
plt.plot(smooth_times, model.simulate(parameters, smooth_times), label='Literature model')
plt.legend()
plt.show()
Explanation: We can also compare the predictions with these values to what we found:
End of explanation
transformation = pints.RectangularBoundariesTransformation(lower, upper)
x0 = [list(parameters) + [0.5]]*3
mcmc = pints.MCMCController(log_posterior, 3, x0,method=pints.HamiltonianMCMC,
                            transformation=transformation)
mcmc.set_max_iterations(100)
chains = mcmc.run()
results = pints.MCMCSummary(
    chains=chains, 
    time=mcmc.time(), 
    parameter_names=["P0", "nu", "k_1", "h", "sigma"])
print(results)
pints.plot.trace(chains, ref_parameters=list(parameters)+[0.75])
plt.show()
Explanation: This is interesting, as this shows that with the sparse sampling points, we cannot retreive the same dynamical behaviour as the literature model. Our model predictions from the MCMC posterior actually give a closer alignment to the data points.
Since this model has sensitivities, we can try Hamiltonian Monte Carlo for inference.
End of explanation
<END_TASK>