In this notebook we set up a simple model and look at sensitivities of states, time-derivatives of states (dot(x)
), and calculated outputs with respect to parameters and initial conditions.
By definition, in Myokit: t0≡0
We can solve this system analytically, to find: x(t)=x0+pqty(t)=2p+(y0−2p)e−t
First, we set up a model and compare the ODE solutions to their analytical equivalent. We check the states (x and y), their time-derivatives (˙x and ˙y), and the outputs (fx and fy).
import matplotlib.pyplot as plt
import myokit
import numpy as np
m = myokit.parse_model('''
[[model]]
e.x = 1.2
e.y = 2.3
[e]
t = 0 bind time
p = 1 / 100
q = 2
r = 3
dot(x) = p * q
dot(y) = 2 * p - y
fx = x^2
fy = r + p^2 + dot(y)
''')
m.validate()
s = myokit.Simulation(m)
d = s.run(8).npview()
t = d.time()
x0 = m.get('e.x').initial_value(True)
y0 = m.get('e.y').initial_value(True)
p = m.get('e.p').eval()
q = m.get('e.q').eval()
r = m.get('e.r').eval()
fig = plt.figure(figsize=(16, 7))
ax = fig.add_subplot(1, 1, 1)
ax.plot(t, d['e.x'], lw=3, label='x')
ax.plot(t, x0 + p * q * t, 'k--')
ax.plot(t, d['e.y'], lw=3, label='y')
ax.plot(t, 2 * p + (y0 - 2 * p) * np.exp(-t), 'k--')
ax.plot(t, d['dot(e.x)'], lw=3, label="x'")
ax.plot(t, np.ones(t.shape) * p * q, 'k--')
ax.plot(t, d['dot(e.y)'], lw=3, label="y'")
ax.plot(t, (2 * p - y0) * np.exp(-t), 'k--')
ax.plot(t, d['e.fx'], lw=3, label='fx')
ax.plot(t, (p * q * t)**2 + 2 * p * q * t * x0 + x0**2, 'k--')
ax.plot(t, d['e.fy'], lw=3, label='fy')
ax.plot(t, r + p**2 + (2 * p - y0) * np.exp(-t), 'k--')
ax.legend()
plt.show()
Sensitivities of x=x0+pqt: ∂x∂p|t=qt∂x∂q|t=pt∂x∂x0|t=1
sens = (['e.x', 'dot(e.x)'], ['e.p', 'e.q', 'init(e.x)'])
s = myokit.Simulation(m, sensitivities=sens)
d, e = s.run(8)
d, e = d.npview(), np.array(e)
t = d.time()
fig = plt.figure(figsize=(16, 7))
ax = fig.add_subplot(1, 2, 1)
ax.plot(t, e[:, 0, 0], lw=3, label='dx/dp')
ax.plot(t, q * t, 'k--')
ax.plot(t, e[:, 0, 2], lw=3, label='dx/dx0')
ax.plot(t, np.ones(t.shape), 'k--')
ax.plot(t, e[:, 1, 0], lw=3, label="dx'/dp")
ax.plot(t, q * np.ones(t.shape), 'k--')
ax.plot(t, e[:, 1, 2], lw=3, label="dx'/dx0")
ax.plot(t, np.zeros(t.shape), 'k--')
ax.legend()
ax = fig.add_subplot(1, 2, 2)
ax.plot(t, e[:, 0, 1], lw=3, color='tab:purple', label='dx/dq')
ax.plot(t, p * t, 'k--')
ax.plot(t, e[:, 1, 1], lw=3, color='tab:brown', label="dx'/dq")
ax.plot(t, p * np.ones(t.shape), 'k--')
ax.legend()
plt.show()
Sensitivities of y=2p+(y0−2p)e−t: ∂y∂p|t=2−2e−t∂y∂y0|t=e−t
sens = (['e.y', 'dot(e.y)'], ['e.p', 'init(e.y)'])
s = myokit.Simulation(m, sensitivities=sens)
d, e = s.run(8)
d, e = d.npview(), np.array(e)
t = d.time()
fig = plt.figure(figsize=(16, 7))
ax = fig.add_subplot(1, 1, 1)
ax.plot(t, e[:, 0, 0], lw=3, label='dy/dp')
ax.plot(t, 2 - 2 * np.exp(-t), 'k--')
ax.plot(t, e[:, 0, 1], lw=3, label='dy/dy0')
ax.plot(t, np.exp(-t), 'k--')
ax.plot(t, e[:, 1, 0], lw=3, label="dy'/dp")
ax.plot(t, 2 * np.exp(-t), 'k--')
ax.plot(t, e[:, 1, 1], lw=3, label="dy'/dy0")
ax.plot(t, -np.exp(-t), 'k--')
ax.legend()
plt.show()
Sensitivities of fx=(pqt)2+2pqtx0+x20: ∂fx∂p|t=2p(qt)2+2qtx0∂fx∂q|t=2p(qt)2+2ptx0∂fx∂x0|t=2pqt+2x0
sens = (['e.fx'], ['e.p', 'e.q', 'init(e.x)'])
s = myokit.Simulation(m, sensitivities=sens)
d, e = s.run(8)
d, e = d.npview(), np.array(e)
t = d.time()
fig = plt.figure(figsize=(16, 7))
ax = fig.add_subplot(1, 2, 1)
ax.plot(t, e[:, 0, 0], lw=3, label='dfx/dp')
ax.plot(t, 2 * p * (q * t)**2 + 2 * q * t * x0, 'k--')
ax.legend()
ax = fig.add_subplot(1, 2, 2)
ax.plot(t, e[:, 0, 1], lw=3, color='tab:orange', label='dfx/dq')
ax.plot(t, 2 * q * (p * t)**2 + 2 * p * t * x0, 'k--')
ax.plot(t, e[:, 0, 2], color='tab:green', lw=3, label='dfx/dx0')
ax.plot(t, 2 * p * q * t + 2 * x0, 'k--')
ax.legend()
plt.show()
Sensitivities of fy=r+p2+(2p−y0)e−t: ∂fy∂p|t=2p+2e−t∂fy∂p|t=1∂fy∂y0|t=−e−t
sens = (['e.fy'], ['e.p', 'e.r', 'init(e.y)'])
s = myokit.Simulation(m, sensitivities=sens)
d, e = s.run(8)
d, e = d.npview(), np.array(e)
t = d.time()
fig = plt.figure(figsize=(16, 7))
ax = fig.add_subplot(1, 1, 1)
ax.plot(t, e[:, 0, 0], lw=3, label='dfy/dp')
ax.plot(t, 2 * p + 2 * np.exp(-t), 'k--')
ax.plot(t, e[:, 0, 1], lw=3, label='dfy/dr')
ax.plot(t, np.ones(t.shape), 'k--')
ax.plot(t, e[:, 0, 2], lw=3, label='dfy/dy0')
ax.plot(t, -np.exp(-t), 'k--')
ax.legend()
plt.show()