# Source code for dynamo.simulation.ODE

```from random import seed, uniform

import anndata
import numpy as np
import pandas as pd

# TODO: import from here in ..estimation.fit_jacobian.py
def hill_inh_func(x, A, K, n, g=0):
Kd = K**n
return A * Kd / (Kd + x**n) - g * x

def hill_inh_grad(x, A, K, n, g=0):
Kd = K**n
return -A * n * Kd * x ** (n - 1) / (Kd + x**n) ** 2 - g

def hill_act_func(x, A, K, n, g=0):
Kd = K**n
return A * x**n / (Kd + x**n) - g * x

def hill_act_grad(x, A, K, n, g=0):
Kd = K**n
return A * n * Kd * x ** (n - 1) / (Kd + x**n) ** 2 - g

[docs]def toggle(ab, t=None, beta=5, gamma=1, n=2):
"""Right hand side (rhs) for toggle ODEs."""
if len(ab.shape) == 2:
a, b = ab[:, 0], ab[:, 1]
res = np.array([beta / (1 + b**n) - a, gamma * (beta / (1 + a**n) - b)]).T
else:
a, b = ab
res = np.array([beta / (1 + b**n) - a, gamma * (beta / (1 + a**n) - b)])

return res

[docs]def Ying_model(x, t=None):
"""network used in the potential landscape paper from Ying, et. al:
https://www.nature.com/articles/s41598-017-15889-2.
This is also the mixture of Gaussian model.
"""
if len(x.shape) == 2:
dx1 = -1 + 9 * x[:, 0] - 2 * pow(x[:, 0], 3) + 9 * x[:, 1] - 2 * pow(x[:, 1], 3)
dx2 = 1 - 11 * x[:, 0] + 2 * pow(x[:, 0], 3) + 11 * x[:, 1] - 2 * pow(x[:, 1], 3)

ret = np.array([dx1, dx2]).T
else:
dx1 = -1 + 9 * x[0] - 2 * pow(x[0], 3) + 9 * x[1] - 2 * pow(x[1], 3)
dx2 = 1 - 11 * x[0] + 2 * pow(x[0], 3) + 11 * x[1] - 2 * pow(x[1], 3)

ret = np.array([dx1, dx2])

return ret

def jacobian_Ying_model(x, t=None):
"""network used in the potential landscape paper from Ying, et. al:
https://www.nature.com/articles/s41598-017-15889-2.
This is also the mixture of Gaussian model.
"""
if len(x.shape) == 2:
df1_dx1 = 9 - 6 * pow(x[:, 0], 2)
df1_dx2 = 9 - 6 * pow(x[:, 1], 2)
df2_dx1 = -11 + 6 * pow(x[:, 0], 2)
df2_dx2 = 11 - 6 * pow(x[:, 1], 2)
else:
df1_dx1 = 9 - 6 * pow(x[0], 2)
df1_dx2 = 9 - 6 * pow(x[1], 2)
df2_dx1 = -11 + 6 * pow(x[0], 2)
df2_dx2 = 11 - 6 * pow(x[1], 2)

J = np.array([[df1_dx1, df1_dx2], [df2_dx1, df2_dx2]])

return J

def hessian_Ying_model(x, t=None):
"""network used in the potential landscape paper from Ying, et. al:
https://www.nature.com/articles/s41598-017-15889-2.
This is also the mixture of Gaussian model.
"""
if len(x.shape) == 2:
H = np.zeros([2, 2, 2, x.shape[0]])
H[0, 0, 0, :] = -12 * x[:, 0]
H[1, 1, 0, :] = -12 * x[:, 1]
H[0, 0, 1, :] = 12 * x[:, 0]
H[1, 1, 1, :] = -12 * x[:, 1]
else:
H = np.zeros([2, 2, 2])
H[0, 0, 0] = -12 * x[0]
H[1, 1, 0] = -12 * x[1]
H[0, 0, 1] = 12 * x[0]
H[1, 1, 1] = -12 * x[1]

return H

def ode_bifur2genes(x: np.ndarray, a=[1, 1], b=[1, 1], S=[0.5, 0.5], K=[0.5, 0.5], m=[4, 4], n=[4, 4], gamma=[1, 1]):
"""The ODEs for the toggle switch motif with self-activation and mutual inhibition (e.g. Gata1-Pu.1)."""

d = x.ndim
x = np.atleast_2d(x)
dx = np.zeros(x.shape)

dx[:, 0] = hill_act_func(x[:, 0], a[0], S[0], m[0], g=gamma[0]) + hill_inh_func(x[:, 1], b[0], K[0], n[0])
dx[:, 1] = hill_act_func(x[:, 1], a[1], S[1], m[1], g=gamma[1]) + hill_inh_func(x[:, 0], b[1], K[1], n[1])

if d == 1:
dx = dx.flatten()

return dx

def jacobian_bifur2genes(
x: np.ndarray, a=[1, 1], b=[1, 1], S=[0.5, 0.5], K=[0.5, 0.5], m=[4, 4], n=[4, 4], gamma=[1, 1]
):
"""The Jacobian of the toggle switch ODE model."""
df1_dx1 = hill_act_grad(x[:, 0], a[0], S[0], m[0], g=gamma[0])
df1_dx2 = hill_inh_grad(x[:, 1], b[0], K[0], n[0])
df2_dx1 = hill_inh_grad(x[:, 0], b[1], K[1], n[1])
df2_dx2 = hill_act_grad(x[:, 1], a[1], S[1], m[1], g=gamma[1])
J = np.array([[df1_dx1, df1_dx2], [df2_dx1, df2_dx2]])
return J

def two_genes_motif_jacobian(x1, x2):
"""This should be equivalent to jacobian_bifur2genes when using default parameters"""
J = np.array(
[
[
0.25 * x1**3 / (0.0625 + x1**4) ** 2 - 1,
-0.25 * x2**3 / (0.0625 + x2**4) ** 2,
],
[
-0.25 * x1**3 / (0.0625 + x1**4) ** 2,
0.25 * x2**3 / (0.0625 + x2**4) ** 2 - 1,
],
]
)
return J

Kd = K**n
return A * n * Kd * x ** (n - 2) * ((n + 1) * x**n - Kd * n + Kd) / (Kd + x**n) ** 3

Kd = K**n
return -A * n * Kd * x ** (n - 2) * ((n + 1) * x**n - Kd * n + Kd) / (Kd + x**n) ** 3

def hessian_bifur2genes(x: np.ndarray, a=[1, 1], b=[1, 1], S=[0.5, 0.5], K=[0.5, 0.5], m=[4, 4], n=[4, 4], t=None):
"""The Hessian of the toggle switch ODE model."""
if len(x.shape) == 2:
H = np.zeros([2, 2, 2, x.shape[0]])
H[0, 0, 0, :] = hill_act_grad2(x[:, 0], a[0], S[0], m[0])
H[1, 1, 0, :] = hill_inh_grad2(x[:, 1], b[0], K[0], n[0])
H[0, 0, 1, :] = hill_inh_grad2(x[:, 0], b[1], K[1], n[1])
H[1, 1, 1, :] = hill_act_grad2(x[:, 1], a[1], S[1], m[1])
else:
H = np.zeros([2, 2, 2])
H[0, 0, 0] = hill_act_grad2(x[0], a[0], S[0], m[0])
H[1, 1, 0] = hill_inh_grad2(x[1], b[0], K[0], n[0])
H[0, 0, 1] = hill_inh_grad2(x[0], b[1], K[1], n[1])
H[1, 1, 1] = hill_act_grad2(x[1], a[1], S[1], m[1])

return H

def ode_osc2genes(x: np.ndarray, a, b, S, K, m, n, gamma):
"""The ODEs for the two gene oscillation based on a predator-prey model."""

d = x.ndim
x = np.atleast_2d(x)
dx = np.zeros(x.shape)

dx[:, 0] = hill_act_func(x[:, 0], a[0], S[0], m[0], g=gamma[0]) + hill_inh_func(x[:, 1], b[0], K[0], n[0])
dx[:, 1] = hill_act_func(x[:, 1], a[1], S[1], m[1], g=gamma[1]) + hill_act_func(x[:, 0], b[1], K[1], n[1])

if d == 1:
dx = dx.flatten()

return dx

def ode_neurogenesis(
x: np.ndarray,
a,
K,
n,
gamma,
):
"""The ODE model for the neurogenesis system that used in benchmarking Monocle 2, Scribe and dynamo (here), original from Xiaojie Qiu, et. al, 2012."""

d = x.ndim
x = np.atleast_2d(x)
dx = np.zeros(shape=x.shape)

dx[:, 0] = (
a[0] * K[0] ** n[0] / (K[0] ** n[0] + x[:, 4] ** n[0] + x[:, 9] ** n[0] + x[:, 11] ** n[0]) - gamma[0] * x[:, 0]
)  # Pax6
dx[:, 1] = (
a[1] * (x[:, 0] ** n[1]) / (K[1] ** n[1] + x[:, 0] ** n[1] + x[:, 5] ** n[1]) - gamma[1] * x[:, 1]
)  # Mash1
dx[:, 2] = a[2] * (x[:, 1] ** n[2]) / (K[2] ** n[2] + x[:, 1] ** n[2]) - gamma[2] * x[:, 2]  # Zic1
dx[:, 3] = (
a[3] * (x[:, 1] ** n[3]) / (K[3] ** n[3] + x[:, 1] ** n[3] + x[:, 7] ** n[3]) - gamma[3] * x[:, 3]
)  # Brn2
dx[:, 4] = (
a[4]
* (x[:, 2] ** n[4] + x[:, 3] ** n[4] + x[:, 10] ** n[4])
/ (K[4] ** n[4] + x[:, 2] ** n[4] + x[:, 3] ** n[4] + x[:, 10] ** n[4])  # Tuj1
- gamma[4] * x[:, 4]
)
dx[:, 5] = (
a[5] * (x[:, 0] ** n[5]) / (K[5] ** n[5] + x[:, 0] ** n[5] + x[:, 1] ** n[5]) - gamma[5] * x[:, 5]
)  # Hes5
dx[:, 6] = a[6] * (x[:, 5] ** n[6]) / (K[6] ** n[6] + x[:, 5] ** n[6] + x[:, 7] ** n[6]) - gamma[6] * x[:, 6]  # Scl
dx[:, 7] = (
a[7] * (x[:, 5] ** n[7]) / (K[7] ** n[7] + x[:, 5] ** n[7] + x[:, 6] ** n[7]) - gamma[7] * x[:, 7]
)  # Olig2
dx[:, 8] = (
a[8] * (x[:, 5] ** n[8] * x[:, 6] ** n[8]) / (K[8] ** n[8] + x[:, 5] ** n[8] * x[:, 6] ** n[8])
- gamma[8] * x[:, 8]  # Stat3
)
dx[:, 9] = a[9] * (x[:, 8] ** n[9]) / (K[9] ** n[9] + x[:, 8] ** n[9]) - gamma[9] * x[:, 9]  # A1dh1L
dx[:, 10] = a[10] * (x[:, 7] ** n[10]) / (K[10] ** n[10] + x[:, 7] ** n[10]) - gamma[10] * x[:, 10]  # Myt1L
dx[:, 11] = a[11] * (x[:, 7] ** n[11]) / (K[11] ** n[11] + x[:, 7] ** n[11]) - gamma[11] * x[:, 11]  # Sox8

if d == 1:
dx = dx.flatten()

return dx

[docs]def neurogenesis(
x,
t=None,
mature_mu=0,
n=4,
k=1,
a=4,
eta=0.25,
eta_m=0.125,
eta_b=0.1,
a_s=2.2,
a_e=6,
mx=10,
):
"""The ODE model for the neurogenesis system that used in benchmarking Monocle 2, Scribe and dynamo (here), original from Xiaojie Qiu, et. al, 2012."""

dx = np.nan * np.ones(shape=x.shape)

if len(x.shape) == 2:
dx[:, 0] = a_s * 1 / (1 + eta**n * (x[:, 4] + x[:, 10] + x[:, 7]) ** n * x[:, 12] ** n) - k * x[:, 0]
dx[:, 1] = a * (x[:, 0] ** n) / (1 + x[:, 0] ** n + x[:, 5] ** n) - k * x[:, 1]
dx[:, 2] = a * (x[:, 1] ** n) / (1 + x[:, 1] ** n) - k * x[:, 2]
dx[:, 3] = a * (x[:, 1] ** n) / (1 + x[:, 1] ** n) - k * x[:, 3]
dx[:, 4] = (
a_e * (x[:, 2] ** n + x[:, 3] ** n + x[:, 9] ** n) / (1 + x[:, 2] ** n + x[:, 3] ** n + x[:, 9] ** n)
- k * x[:, 4]
)
dx[:, 5] = a * (x[:, 0] ** n) / (1 + x[:, 0] ** n + x[:, 1] ** n) - k * x[:, 5]
dx[:, 6] = a_e * (eta**n * x[:, 5] ** n) / (1 + eta**n * x[:, 5] ** n + x[:, 7] ** n) - k * x[:, 6]
dx[:, 7] = a_e * (eta**n * x[:, 5] ** n) / (1 + x[:, 6] ** n + eta**n * x[:, 5] ** n) - k * x[:, 7]
dx[:, 8] = (
a * (eta**n * x[:, 5] ** n * x[:, 6] ** n) / (1 + eta**n * x[:, 5] ** n * x[:, 6] ** n) - k * x[:, 8]
)
dx[:, 9] = a * (x[:, 7] ** n) / (1 + x[:, 7] ** n) - k * x[:, 9]
dx[:, 10] = a_e * (x[:, 8] ** n) / (1 + x[:, 8] ** n) - k * x[:, 10]
dx[:, 11] = a * (eta_m**n * x[:, 7] ** n) / (1 + eta_m**n * x[:, 7] ** n) - k * x[:, 11]
dx[:, 12] = mature_mu * (1 - x[:, 12] / mx)
else:
dx[0] = a_s * 1 / (1 + eta**n * (x[4] + x[10] + x[7]) ** n * x[12] ** n) - k * x[0]
dx[1] = a * (x[0] ** n) / (1 + x[0] ** n + x[5] ** n) - k * x[1]
dx[2] = a * (x[1] ** n) / (1 + x[1] ** n) - k * x[2]
dx[3] = a * (x[1] ** n) / (1 + x[1] ** n) - k * x[3]
dx[4] = a_e * (x[2] ** n + x[3] ** n + x[9] ** n) / (1 + x[2] ** n + x[3] ** n + x[9] ** n) - k * x[4]
dx[5] = a * (x[0] ** n) / (1 + x[0] ** n + x[1] ** n) - k * x[5]
dx[6] = a_e * (eta**n * x[5] ** n) / (1 + eta**n * x[5] ** n + x[7] ** n) - k * x[6]
dx[7] = a_e * (eta**n * x[5] ** n) / (1 + x[6] ** n + eta**n * x[5] ** n) - k * x[7]
dx[8] = a * (eta**n * x[5] ** n * x[6] ** n) / (1 + eta**n * x[5] ** n * x[6] ** n) - k * x[8]
dx[9] = a * (x[7] ** n) / (1 + x[7] ** n) - k * x[9]
dx[10] = a_e * (x[8] ** n) / (1 + x[8] ** n) - k * x[10]
dx[11] = a * (eta_m**n * x[7] ** n) / (1 + eta_m**n * x[7] ** n) - k * x[11]
dx[12] = mature_mu * (1 - x[12] / mx)

return dx

def hsc():
pass

[docs]def state_space_sampler(ode, dim, seed_num=19491001, clip=True, min_val=0, max_val=4, N=10000):
"""Sample N points from the dim dimension gene expression space while restricting the values to be between min_val and max_val. Velocity vector at the sampled points will be calculated according to ode function."""

seed(seed)
X = np.array([[uniform(min_val, max_val) for _ in range(dim)] for _ in range(N)])
Y = np.clip(X + ode(X), a_min=min_val, a_max=None) if clip else X + ode(X)

return X, Y

[docs]def Simulator(motif="neurogenesis", seed_num=19491001, clip=True, cell_num=5000):
"""Simulate the gene expression dynamics via deterministic ODE model

Parameters
----------
motif: str (default: "neurogenesis")
Name of the network motif that will be used in the simulation.
clip: bool (default: True)
Whether to clip data points that are negative.
cell_num: int (default: 5000)
Number of cells to simulate.

Returns
-------
an Annodata object containing the simulated data.
"""

if motif == "toggle":
X, Y = state_space_sampler(
ode=toggle,
dim=2,
seed_num=seed_num,
min_val=0,
max_val=6,
N=cell_num,
clip=clip,
)
gene_name = np.array(["X", "Y"])
elif motif == "neurogenesis":
X, Y = state_space_sampler(
ode=neurogenesis,
dim=13,
seed_num=seed_num,
min_val=0,
max_val=6,
N=cell_num,
clip=clip,
)

gene_name = np.array(
[
"Pax6",  #
"Mash1",  #
"Brn2",
"Zic1",
"Tuj1",
"Hes5",  #
"Scl",  #
"Olig2",  #
"Stat3",
"Myt1L",
"Alhd1L",
"Sox8",
"Maturation",
]
)
elif motif == "twogenes":
X, Y = state_space_sampler(ode=ode_bifur2genes, dim=2, min_val=0, max_val=4, N=cell_num, clip=clip)
gene_name = np.array(["Pu.1", "Gata.1"])
elif motif == "Ying":
X, Y = state_space_sampler(ode=Ying_model, dim=2, min_val=-3, max_val=3, N=cell_num, clip=clip)
gene_name = np.array(["X", "Y"])

var = pd.DataFrame({"gene_short_name": gene_name})  # use the real name in simulation?
var.set_index("gene_short_name", inplace=True)

# provide more annotation for cells next:
cell_ids = ["cell_%d" % (i) for i in range(cell_num)]  # first n_traj and then steps
obs = pd.DataFrame({"Cell_name": cell_ids})
obs.set_index("Cell_name", inplace=True)

layers = {"velocity": Y - X}  # ambiguous is required for velocyto

adata = anndata.AnnData(X.copy(), obs.copy(), var.copy(), layers=layers.copy())

# remove cells that has no expression