import numpy as np
from random import uniform, seed
import anndata
import pandas as pd
[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"""
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
[docs]def two_genes_motif(x, t=None, a1=1, a2=1, b1=1, b2=1, k1=1, k2=1, S=0.5, n=4):
"""The ODE model for the famous Pu.1-Gata.1 like network motif with self-activation and mutual inhibition."""
dx = np.nan * np.ones(x.shape)
if len(x.shape) == 2:
dx[:, 0] = a1 * x[:, 0] ** n / (S ** n + x[:, 0] ** n) + b1 * S ** n / (S ** n + x[:, 1] ** n) - k1 * x[:, 0]
dx[:, 1] = a2 * x[:, 1] ** n / (S ** n + x[:, 1] ** n) + b2 * S ** n / (S ** n + x[:, 0] ** n) - k2 * x[:, 1]
else:
dx[0] = a1 * x[0] ** n / (S ** n + x[0] ** n) + b1 * S ** n / (S ** n + x[1] ** n) - k1 * x[0]
dx[1] = a2 * x[1] ** n / (S ** n + x[1] ** n) + b2 * S ** n / (S ** n + x[0] ** n) - k2 * x[1]
return dx
def two_genes_motif_jacobian(x1, x2):
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
[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, 2011."""
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
[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
-------
adata: :class:`~anndata.AnnData`
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,
)
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,
)
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=two_genes_motif, dim=2, min_val=0, max_val=4, N=cell_num)
gene_name = np.array(["Pu.1", "Gata.1"])
elif motif == "Ying":
X, Y = state_space_sampler(ode=Ying_model, dim=2, clip=clip, min_val=-3, max_val=3, N=cell_num)
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
adata = adata[adata.X.sum(1) > 0, :] if clip else adata
return adata