Bayesian Inference Algorithm for Estimating Heterogeneity of Regulatory Mechanisms Based on Single-Cell Data

Abstract

The rapid progress in biological experimental technologies has generated a huge amount of experimental data to investigate complex regulatory mechanisms. Various mathematical models have been proposed to simulate the dynamic properties of molecular processes using the experimental data. However, it is still difficult to estimate unknown parameters in mathematical models for the dynamics in different cells due to the high demand for computing power. In this work, we propose a population statistical inference algorithm to improve the computing efficiency. In the first step, this algorithm clusters single cells into a number of groups based on the distances between each pair of cells. In each cluster, we then infer the parameters of the mathematical model for the first cell. We propose an adaptive approach that uses the inferred parameter values of the first cell to formulate the prior distribution and acceptance criteria of the following cells. Three regulatory network models were used to examine the efficiency and effectiveness of the designed algorithm. The computational results show that the new method reduces the computational time significantly and provides an effective algorithm to infer the parameters of regulatory networks in a large number of cells.

1. Introduction

The fast progress in biological experimental technologies in recent years has generated a huge amount of experimental data to investigate the molecular regulatory mechanisms inside the cell [1,2,3,4]. Among them, single-cell technologies have been used to quantify the expression profiles and protein activities in different single cells at the same time [5,6]. These quantitative and qualitative datasets provide rich information for studying the heterogeneity of regulatory mechanisms in different cells and in different patients for disease therapy. However, there are substantial challenges in illustrating the origin of heterogeneity and to describe the propagation of heterogeneity in cellular processes [7,8].

Mathematical modeling is a powerful method to investigate the diverse dynamic observations in a large number of cells. For time-lapse data, the two-stage modeling method uses a single mathematical model with different model parameter values to simulate the diverse dynamics in single cells [9]. The heterogeneity of biological processes is demonstrated by the distributions of the model parameters, which may be overestimated. To address this issue, the global two-stage method uses a local measure for the distributions of the model parameters [10]. Another widely used method is the mixed-effect model, which derives the likelihood functions of the model parameters by employing the first-order conditional estimation [11,12] or the stochastic approximation expectation–maximization algorithm [11]. In addition, the maximum entropy method has been used to analyze the single-cell snapshot data [13]. Furthermore, the likelihood function has been used to identify the variability in biological processes and the distributions of the model parameters [14,15]. More related research studies can be found in the review articles published in recent years [16,17].

In this study, we discuss inference methods for estimating the unknown model parameters based on single-cell time-lapse data. Thus, we estimate a number of mathematical models using the same model equations, but with different sets of model parameters. For the inference of the unknown model parameters, the Bayesian statistical methods and optimization algorithms are two major types of widely used methods [18,19,20]. In recent years, machine learning methods have been used to estimate unknown model parameters [21]. Among these methods, the Bayesian inference methods are able to obtain the distributions of the model parameters and also have the ability to analyze noisy datasets with better accuracy. However, the classic Markov chain Monte Carlo (MCMC) method needs the likelihood function to assess the accuracy of the model parameters, which are difficult to apply to complex systems without explicit likelihood functions. To overcome this limitation, the approximate Bayesian computation (ABC) algorithms have been proposed to measure the accuracy of the model parameters by comparing the model simulations with experimental data directly [22]. The first proposed ABC was the ABC rejection algorithm, which needs a large computing time because this method has no learning step. More effective algorithms have been designed to speed up the convergence rate, for example ABC-MCMC and sequential Monte Carlo (SMC) ABC [23,24,25]. In recent year, more measures have been designed to assess the accuracy of the model simulations [26]; the early rejection algorithms have been developed to reduce the computing time [27].

The ABC-SMC algorithm uses the adaptive transitional kernel functions to accelerate the acceptance rate [24]. However, this may be difficult for the model solutions close to the experimental data if the data are noisy. Thus, the selection of the tolerance threshold values is critical for the successful applications of ABC-SMC. One option for selecting the threshold values is to use the simulation errors of the accepted samples to construct the adaptive tolerance threshold [23,28]. However, it is not easy to use these adaptive approaches to estimate the model parameters using noisy data. For the application problems, a simple approach is to choose the threshold values by manual adjustment. Since the single-cell time-lapse datasets may include observations from a large number of single cells, it is not practical to adjust the threshold values for each single cell.

We conducted initial computations to infer the parameters in a population of models [29]. This study used the existing ABC-SMC to infer the parameters of each model, which requires a huge computing time. To address the identified challenges in our initial study, we designed a new algorithm in this work to reduce the computational time of ABC-SMC. The innovation of this study is dividing the dataset with many single cells into a few clusters. Each cluster has a number of single cells whose expression profiles are close to each other. After obtaining the estimates of the parameters for the first cell in each cluster, we designed an adaptive method to use these estimates to construct the threshold values for the other cells in the same cluster. Three test system models were used to evaluate the efficiency of the proposed algorithm.

2. Materials and Methods

2.1. Mathematical Model

In this work, we studied the inference methods for estimating the unknown parameters of the following system:

$$\begin{array}{ccc}\hfill \frac{d{x}_i}{dt}& =& f_i(x_1,x_2,\cdots ,x_n,B_1,B_2,\dots ,B_m),i=1,2,\dots ,n,\hfill \end{array}$$

(1)

where ${(x_1\left(t\right),\dots ,x_n\left(t\right))}^{\top }$ are the system states at time t, $(B_1,\dots ,B_m)$ are model parameters, and functions $(f_1,\dots ,f_n)$ describe the evolutions of the system state. If the value of $B_i$ ($i=1,\dots ,m$) is constant, System (1) is the traditional ordinary differential equations (ODEs). Here, we studied the case that $B_i$ are random variables following a joint distribution. Note that, unlike stochastic differential equations, the value of $B_i$ is not a sample of a stochastic process over time. To simulate the dynamics in single cells, it was assumed that the model for each single cell has a particular value of parameter $B_i$, and the model of different single cells has different values. These values are samples of the particular joint distribution.

In the inference study, we first generated a sample of parameters $(b_1,\dots ,b_m)$ from a prior distribution $\pi \left(\theta \right)$. Since we did not have any information about the prior distribution, it was assumed that this distribution is a uniform distribution over the interval $[a,b]$. The values of a and b are selected in the inference steps. In addition, it was assumed that random variables $B_i$ are independent of each other. Once a sample $(b_1,\dots ,b_m)$ is obtained, we used the differential equation system (1), in which random variable $B_i$ is replaced by sample $b_i$, to simulate the molecular dynamics in single cells.

2.2. Approximation Bayesian Computation

This subsection briefly introduces the widely used ABC-SMC algorithm, which was used in this work to infer the model parameters. We also discuss the issues in the applications of this algorithm, which is given below [30]. More detailed information regarding the implementations of this algorithm can be found in [22].

The discrepancy between the simulation and experimental data is the mean-squared relative error, calculated by

$$\rho (X,Y)=\frac{1}{mT}\sum _{i=1}^m\sum _{j=0}^T{\left(\frac{x_{ij}-y_{ij}}{x_{ij}}\right)}^2,$$

(2)

where $y_{ij}$ and $x_{ij}$ are the experimental data and model simulation of the i-th variable $x_i$ at time point $t_j$, respectively.

One question in the application of ABC-SMC is the choice of the proper threshold values ${ϵ}_1,\dots ,{ϵ}_K$. Although we may be able to manually select these values for the experimental data in each single cell, it is time consuming to choose these values manually when the cell number is large. Thus, we needed to design a technique to select the threshold values effectively for a large number of cells.

2.3. ABC-SMC with Adaptive Tolerance Threshold

To improve the efficiency, a natural idea is to use the estimated parameters for the first cell to construct the threshold values and prior distribution for the cell that is close to the first cell. We can use particles in the last generation (i.e., the K-th generation in Algorithm 1) to develop the prior distribution and use the corresponding discrepancies in Step 3.4 in Algorithm 1 to design the tolerance threshold.

Algorithm 1 ABC-SMC algorithm.

1: Input information: experimental data $X=\{x_0,x_1,\dots ,x_T\}$ at time point $t=[t_0,t_1,\dots ,t_T]$; prior distribution of parameters $\pi \left(\theta \right)$, tolerance threshold values ${ϵ}_1,\dots ,{ϵ}_K$, where K is the number of generations. 2: The first generation $k=1$ For particle $j=1,\dots ,M$    (1) Use the prior distribution to generate a sample ${\theta }^{*}\sim \pi \left(\theta \right)$.    (2) Simulate the model with parameters ${\theta }^{*}$ to obtain the simulation data $X^{*}\sim$ Model(${\theta }^{*}$).    (3) Calculate the discrepancy between simulation data and experimental data $\rho (X,X^{*})$.    (4) If $\rho (X,X^{*})>{ϵ}_1$, reject the sample. Otherwise, accept this sample as ${\theta }_{j,1}={\theta }^{*}$.    (5) Set the same weight to each particle as $w_{j,1}=\frac{1}{M}$. Set the variance of the particles in the first generation ${\sigma }_k^2=2Var\left({\theta }_{1:M,1}\right)$. 3: For generations $k=2,\dots ,K$ For particle $j=1,\dots ,M$:    (1) Based on the accepted samples in the previous generation, select a sample ${\theta }^{*}$ using the weights of the previous generation $w_{1:M,k-1}$.    (2) Generate a new sample ${\theta }^{**}\sim q\left(\theta \right|{\theta }^{*},{\sigma }_{k-1})$ using the proposed distribution.    (3) Simulate the model with parameter ${\theta }^{**}$ to obtain the simulation data $X^{**}\sim$ Model(${\theta }^{**}$).    (4) Calculate the discrepancy between the simulation data and experimental data $\rho (X,X^{**})$.    (5) If $\rho (X,X^{**})>{ϵ}_k$, reject the sample. Otherwise, accept this sample as ${\theta }_{j,k}={\theta }^{**}$. Calculate the new weights for the accepted samples in generation k:    (1) The new weight is $w_{j,k}={\displaystyle \frac{\pi \left({\theta }_{j,k}\right)}{{\sum }_{i=1}^M{w}_{i,k-1}q\left({\theta }_{i,k}\right|{\theta }_{j,k},{\sigma }_{k-1})}}$. The new variance is ${\sigma }_k^2=2Var\left({\theta }_{1:M,k}\right)$. 4: Output: the accepted samples in the K-th generation ${\theta }_{1:M,K}$ and the corresponding discrepancies ${\rho }_{1:M,K}$.

However, this idea cannot be applied directly if the experimental data in two cells have a large distance. For example, these two cells may belong to different cell types, or they may be at different developmental stages. Mathematically, the observed dynamics may show different trends. In addition, after obtaining the estimated parameters for a number of cells, we need to select an optimal estimate to match the dynamics of the next cell.

To address these issues, two improvements of Algorithm 1 are proposed in this work. First, we divided all observed cells into a number of clusters based on the distance between each pair of cells. Since we needed to develop the initial information (i.e., the prior distribution and threshold values) for the first cell in each cluster, the cluster number should not be very small (i.e., the difference between cells in one cluster may be large) or very large (i.e., computational inefficiency). In this study, the cluster number for the three test systems was selected in the range between 10 and 15, depending on the variations of the observations in the dataset.

We first used the mean values of all cells in each cluster to determine the centroid of that cluster and then found the first cell that had the smallest distance to the cluster centroid. Alternatively, we can use medoid clustering for grouping cells and determining the first cell. In medoid clustering, we chose an actual data point as the centroid of a cluster. This is different from the k-means clustering, which uses the mean values as the centroid of a cluster. Based on the first cell, we ranked all the cells in this cluster based on their distances to the first cell. We first used Algorithm 1 to infer the unknown parameters of the first cell, which were then employed to determine the initial information of the second cell. For the k-th cell ($k>2$), we found a cell from the cells $1,2,\dots ,k-1$ that has the smallest distance to the k-th cell and, then, used the estimated parameters of that cell to determine the initial information of the k-th cell.

Except for the first cell, the prior distribution of each parameter ${\theta }_i$ was assumed to follow the uniform distribution $\pi \left(\theta \right)\sim U(W_{min},W_{max})$. Assume that the estimated values of this parameter in the previous cell are ${\theta }_{K,i}$ in the final K-th generation. The boundaries of the uniform distribution are

$$W_{min}=k_{1}min\left\{{\theta }_{1:M,K}\right\},W_{max}=k_{1}max\left\{{\theta }_{1:M,K}\right\},$$

(3)

where parameter $k_1$ is associated with the distance between these two cells. For the tolerance threshold, using the M discrepancy values ${\rho }_{1:M,K}$ in the final generation of the previous cell, we set the first tolerance threshold as

$${ϵ}_1=k_{2}max\left\{{\rho }_{1:M,K}\right\},$$

(4)

where parameter $k_2$ is associated with the distance between these two cells. For the following generations, we used the discrepancy values of the previous generation to determine the tolerance threshold adaptively. The major steps of the proposed method are given in Algorithm 2.

Algorithm 2 ABC-PMC with adaptive tolerance threshold (ABC-CPMC).

1: Initiation: observation data $X=\{X_1,X_2,\dots ,X_N\}$ of N cells. 2: Cluster these N cells into C clusters, and find the centroid of each cluster. Each cluster has $N_i$ cells ($i=1,\dots ,C$). 3: For clusters $i=1,\dots ,C$: (1) Based on the centroid $X_i^{*}$ of this cluster, find the first cell of this cluster by ${min}_j\rho (X_{ij}-X_i^{*})$ ($j=1,\dots ,N_i$). Denote the data of the first cell as $X_{i1}$. (2) Rank all $N_i$ cells in the i-th cluster as $X_{i1},X_{i2},\dots ,X_{i{N}_i}$ based on their distances from the first cell. (3) For the first cell $X_{i1}$, use Algorithm 1 and prior ${\theta }^{*}\sim \pi \left(\theta \right)$ to infer the model parameters. Denote the inferred parameter as ${\theta }_{1:M,K}^{i1}$. (4) For the second cell $X_{i2}$, construct the prior as ${\theta }^{*}\sim U(a,b)$ using Formula (3). The tolerance threshold ${ϵ}_1$ is determined by using Formula (4). (5) Use Algorithm 1 to estimate the parameters ${\theta }_{1:M,K}^{i2}$ of the second cell. (6) For cells $k=3,\dots ,N_i$: (a) Find a cell from the cells $1,2,\dots ,k-1$ that has the smallest distance to the k-th cell in these $k-1$ cells. Denote the particles in the final generation of this cell as ${\theta }_{1:M,K}^{i(k-1)*}$. (b) Construct the prior as ${\theta }^{*}\sim U(a,b)$ using Formula (3). The tolerance threshold ${ϵ}_1$ is determined by using Formula (4). (c) Use Algorithm 1 to estimate the parameters ${\theta }_{1:M,K}^{ik}$ of the k-th cell.

Note that we can also use the proposed technique in Algorithm 2 for ranking cells in a cluster in Steps 3.1 and 3.2 to rank all clusters based on the centroid of each cluster. In this way, the prior distribution of the first cell in a cluster may be obtained by the estimated parameters of the neighboring cluster.

3. Results and Discussion

3.1. Gene Network with One Steady State

We first considered a mathematical model for studying the dynamics of a single gene [31]. This gene has a positive regulation of itself, and this regulation is described by a Michaelis–Menten function. The produced mRNA $\left(x\right)$ from the first equation below will synthesize the proteins $\left(y\right)$ in the second equation. The model is given by

$$\begin{array}{cc}\hfill \frac{dx}{dt}& =\frac{ay}{1+y}-k_{1}x,\hfill \\ \hfill \frac{dy}{dt}& =bx-k_{2}y,\hfill \end{array}$$

(5)

where a and $k_1$ are the synthesis rate and degradation rate of the mRNA and b and $k_2$ are the synthesis rate and degradation rate of the protein, respectively. These four rates are assumed to follow the Gaussian distributions, namely a and $b\sim N(2,1)$, as well as $k_1$ and $k_2\sim N(1,{0.5}^2)$. We used a fixed initial condition: $(x_0,y_0)=(2,3)$.

We tested this system by using 500 simulations [29]. In this study, we extended the simulation number to 1000, shown in Figure 1. Since a Gaussian distribution may generate samples with a wide range of values, we restricted the samples values of a and b in the interval $(1,3)$ and those of $k_1$ and $k_2$ in the interval $(0.5,1.5)$. In addition, simulations with large outliers (i.e., $x>3.5$ or $y>7$) were removed.

We first divided the 1000 generated simulations in Figure 1 into 15 clusters using the command kmeans.m in MATLAB. The number of cluster was selected to avoid mixing the cells in different groups. The minimal distance from cells outside a cluster to the centroid of that cluster is larger than the maximal distance from cells inside the cluster to their centroid. For each simulation, we averaged the values of x and y over the 21 observation time points (i.e., $t=0,1,\dots ,20$) and, then, drew the scatter plot of these 1000 simulations in Figure 2A. This shows that the simulations of the same cluster have small distances to each other. All the clusters are connected, and Clusters 3 and 5 are located at the two ends of the spectrum. Note that this figure gives a two-dimensional representation of the 21-dimensional data. The mixture of the points in the figure does not mean these points mix in the two-dimensional space.

When inferring the unknown parameters, the prior distribution was assumed to be a uniform distribution $\pi \left(\theta \right)\sim U(W_{min},W_{max})$ for each parameter for the first cell in each cluster. The values of $W_{min}$ and $W_{max}$ are $(1,1,0.5,0.5)$ and $(3,3,1.5,1.5)$ for $(a,b,k_1,k_2)$, respectively. The tolerance threshold values ${ϵ}_i$ ($i=1,\dots ,10$) were selected manually as small as possible for the first cell.

Figure 3 provides the simulated frequencies of the estimated parameters and those of the exact parameters of Model (5). The simulated distributions have very good accuracy with respect to the exact ones. However, the variances of the inferred distributions are larger than those of the exact ones because the prior distributions were assumed to be the uniform distributions. The simulated frequencies can be used as the posterior distributions for further analysis.

3.2. Genetic Toggle Switch Showing Bistability

The second test system is the genetic toggle switch model for the expression of the two genes [31,32]. The repressor genes $\lambda$CI (u) and LacR (v) negatively regulate the other gene. The mathematical model is given by

$$\begin{array}{cc}\hfill \frac{du}{dt}& ={\alpha }_1+\frac{{\beta }_1}{1+v^3}-\left(1+\frac{s}{1+s}\right)u,\hfill \\ \hfill \frac{dv}{dt}& ={\alpha }_2+\frac{{\beta }_2}{1+u^3}-v,\hfill \end{array}$$

(6)

where ${\alpha }_1$ and ${\alpha }_2$ are the rate constants of the basal synthesis of the two genes and ${\beta }_1$ and ${\beta }_2$ are the rate constants to realize the negative regulations from the other gene, and s is used to realize genetic switching regulated by protein RecA.

We tested this model (6) by using 500 simulations [29]. In this study, we extended the simulation number to 1000, shown in Figure 1. The fixed initial condition was $(u_0,v_0)=(4.1341,0.2558)$. The system parameters were assumed to obey Gaussian distributions, namely ${\alpha }_1,{\alpha }_2\sim N(0.2,{0.01}^2)$, ${\beta }_1,{\beta }_2\sim N(4,{0.5}^2)$, and $s\sim N(2.135,{0.1}^2)$. This model shows two steady states, but each simulation has only one of these two stable states. The distributions of the steady states at $t=50$ obtained by different model parameters are given in Figure 4C,D.

We first divided the 1000 generated simulations in Figure 1 into 12 clusters using kmeans.m in MATLAB. Similar to the first test system, the number of cluster was selected to avoid mixing the cells in different groups. For each simulation, we averaged the values of u and v over the 101 observation time points (i.e., $t=0,1,2,\dots ,100$) and, then, drew the scatter plot of these 1000 simulations in Figure 2B. The two separate groups of clusters are clearly shown, which indicates the two stable steady states of the system. This result suggests the importance of clustering. When a population of cells has different cell types, the cellular dynamics may be distinct in different cell types. Figure 2B clearly gives the characterization of the bistability property of the gene network. We can use different initial information to estimate the unknown parameters for the cells in these two groups.

Based on the clusters shown in Figure 2B, we next inferred the model parameters for each cell using our proposed algorithm. For the first cell of each cluster, the tolerance threshold values of 10 generations were set as $ϵ=\{1,0.5,0.1,0.05,0.01,0.005,0.002$, $0.001,0.0005,0.0002\}$ initially.

Figure 5 provides the simulated frequencies of estimated parameters and those of the exact parameters of Model (6). Figure 5B–D show that the simulated distributions have very good accuracy with respect to the exact ones for ${\beta }_1$, ${\alpha }_2$ and ${\beta }_2$. However, Figure 5A suggests that the accuracy of the estimates for parameter s is not as good as those of the other three parameters. The reason for this low accuracy may be that parameter s is not zero only in a short time period. The simulated frequencies can be used as the posterior distributions for further analysis.

3.3. MAP Kinase Pathway for Efficiency Test

We tested the accuracy of our proposed algorithm by using two small-scale models of genetic regulations. The key question is whether the new method can achieve substantial improvement of the efficiency over the existing methods. To answer this question, we next considered the mathematical model of the MAP kinase pathway, which is one of the most prominent components of the pathways that control cell proliferation, differentiation, and apoptosis. Based on the observation data in single cells, this work studied a subnetwork of the MAP kinase pathway using the activated Raf protein as the input [14]. The activated Raf protein, denoted as pRaf, activates MEK proteins by phosphorylation, leading to double-phosphorylated MEK (ppMEK). The activated MEK protein then activates the ERK proteins by phosphorylation, leading to single-phosphorylated ERK (pERK) and double-phosphorylated ERK (ppERK). Meanwhile, the phosphatases MEK-P’ase (MEKP) and ERK-P’ase (ERKP) can deactivate ppMEK and ppERK, respectively [33,34]. The detailed process of the kinase activations consists of seven sets of biochemical reactions [35], which are given below. The detailed model of the ODEs is provided in the Supplementary Materials.

$R_1$: pRaf + MEK $\underset{d_1}{\overset{a_1}{\rightleftharpoons }}$ pRaf − MEK $\stackrel{k_1}{\to }$ pRaf + ppMEK,

$R_2$: ppMEK + MEKP $\underset{d_2}{\overset{a_2}{\rightleftharpoons }}$ ppMEK − MEKP $\stackrel{k_2}{\to }$ pMEK + MEKP,

$R_3$: pMEK + MEKP $\underset{d_3}{\overset{a_3}{\rightleftharpoons }}$ pMEK − MEKP $\stackrel{k_3}{\to }$ MEK + MEKP,

$R_4$: ppMEK + ERK $\underset{d_4}{\overset{a_4}{\rightleftharpoons }}$ ppMEK − ERK $\stackrel{k_4}{\to }$ ppMEK + pERK,

$R_5$: ppMEK + pERK $\underset{d_5}{\overset{a_5}{\rightleftharpoons }}$ ppMEK − pERK $\stackrel{k_5}{\to }$ ppMEK + ppERK,

$R_6$: ppERK + ERKP $\underset{d_6}{\overset{a_6}{\rightleftharpoons }}$ ppERK − ERKP $\stackrel{k_6}{\to }$ pERK + ERKP,

$R_7$: pERK + ERKP $\underset{d_7}{\overset{a_7}{\rightleftharpoons }}$ pERK − ERKP $\stackrel{k_7}{\to }$ ERK + ERKP.

The initial protein concentrations were [MEK] = 1.4, [ERK] = 0.96, [MEKP] = 0.7, and [ERKP] = 0.48 [36,37], and the concentrations of the other proteins were zero. As the signal input, the value of pRaf transmitted into the ODE solver in MATLAB used a linear interpolation based on the experimental data used in [36]. The MEK activities and ERK activities were measured from Figure 2 in a recent single-cell study [14]. The experimental data in [14] provide measured values only at seven time points in $[0,20]$. In addition, the measured values are the number of protein molecules, rather than the concentration, in the previous studies. To consolidate different experimental data, we used the kinase activities at 5 min to normalize the data at other time points. Cubic spline interpolation was used to estimate the missing experimental data by using the measured kinase activities.

In the first step, we estimated one set of the model parameters using the average MEK and ERK activities from the single cells [14]. We used ABC-PMC (Algorithm 1) to infer the 21 model parameters. Figure 6 gives the average MER and ERK activities and the simulation using the estimated model parameters. The inferred network dynamics is very consistent with the experimental data.

Table 1. Estimated model parameters for the MAP kinase pathway. Estimated value: the estimated value for each parameter. STD: standard deviation of the 100 estimates in the final generation. $[W_{min},W_{max}]$: prior distribution of each parameter.

Kinetic rates	${\mathit{a}}_1$	${\mathit{a}}_2$	${\mathit{a}}_3$	${\mathit{a}}_4$	${\mathit{a}}_5$	${\mathit{a}}_6$	${\mathit{a}}_7$
Estimated value	66.0452	0.9584	0.0121	15.3943	35.7607	5.7297	5.0556
STD	13.7871	0.6568	0.0031	0.3849	0.9216	0.1316	2.0026
$W_{min}$	20	0	0	10	25	1	0
$W_{max}$	120	5	0.05	20	50	10	15
Kinetic rates	$d_1$	$d_2$	$d_3$	$d_4$	$d_5$	$d_6$	$d_7$
Estimated value	0.1176	60.4114	35.3809	36.2956	16.3216	44.2443	2.4649
STD	0.08	3.2385	1.2341	2.0504	0.6656	1.5363	0.0958
$W_{min}$	0	30	25	25	10	30	1
$W_{max}$	0.5	80	50	50	25	60	5
Kinetic rates	$k_1$	$k_2$	$k_3$	$k_4$	$k_5$	$k_6$	$k_7$
Estimated value	25.248	12.1591	5.3689	59.3748	29.3347	28.6955	27.5407
STD	0.7066	0.4916	0.4452	3.551	1.6659	1.2889	0.553
$W_{min}$	20	5	1	40	20	20	20
$W_{max}$	50	20	10	80	50	50	50

gives the inferred parameter values and corresponding standard deviation based on the 100 particles of the last generation in Algorithm 1.

Based on the inferred model parameters, we next generated simulations using the perturbed model parameters. It was assumed that each model parameter follows a Gaussian distribution whose mean and standard deviation are the values shown in Supplementary Table S1. After generating a sample of the model parameters, we obtained a simulation of the system and, then, examined whether the simulation was within the observed range of kinase activities [14]. If a simulation was well beyond the observation range, it was discarded. In this way, we generated 1000 simulations, which were treated as the observation data for inferring the model parameters in the following step.

Since the model has 15 variables, it is not easy to reduce the dimension of the system. We used Ward’s linkage method [38] to conduct hierarchical clustering of these 1000 simulations. All simulations were divided into 12 clusters. Figure 7 gives the simulated ppMEK and ppERK activities in three clusters. It shows that the differences of the kinase activities in each cluster are small. However, the variations between different clusters are not small.

Using the simulations in the 12 clusters, we used the proposed algorithm to infer the distributions of the model parameters in the final step. To obtain a reasonable prior distribution for each parameter, we studied the sensitivity of that parameter by using the derivatives of the solution with respect to it. We used the package using the iterative approximations based on the directional derivatives [39]. The sensitivity analysis results given in Supplementary Figure S1 suggest that a part of the parameters (e.g., $k_3$) had a weak influence on the system dynamics. Note that the derivatives of ppMEK and ppERK are almost negative with respect to a small neighborhood of parameters such as $d_1=0.1176$. The negative perturbation reaches the maximum, which suggests that we can reduce the simulation when it is larger by adjusting the value of $d_1$. On the other hand, it is convenient to adjust the parameters such that the model solutions do not exceed the reasonable range at each time point. It was assumed that the prior distribution of the parameters follows $\pi \left(\theta \right)\sim U(W_{min},W_{max})$. The values of $W_{min}$ and $W_{max}$ are shown in Table 1.

Figure 8 gives the distributions for two sets of inferred model parameters $(a_i,d_i,k_i)$. Compared with the generated model parameters, the inferred parameters have a relatively larger range of values. The means of the inferred parameters are consistent with those of the generated parameters. However, the variance of the inferred parameters is larger than those of the true parameters.

Finally, we calculated the speedup of the proposed algorithm, which is the key advantage of the new method. We used five clusters to measure the computational time to infer the model parameters. The speedup is defined as the ratio of the computational time of ABC-PMC (Algorithm 1) to that of our proposed algorithm (Algorithm 2). Figure 9 shows that the speedup for these five clusters is around 100, which suggests that our new algorithm achieved a substantial improvement over the existing inference methods.

4. Conclusions

This study designed a novel method for estimating the unknown parameters in a large number of models based on the experimental single-cell data. We clustered single cells into a number of groups based on the distances between each pair of cells. We used simulations to select the initial information (i.e., prior distribution and tolerance threshold) to estimate the unknown parameter for the model of the first cell in each cluster. The inference results from the first step were used to develop the initial information for the following cell. This method was repeated to estimate the parameters of the other cells in one cluster. Three network models were used to evaluate the efficiency and accuracy of this new algorithm. In particular, the genetic toggle switch model was used to show the function of clustering, and the MAP kinase pathway model was employed to demonstrate the computing efficiency. The inference results of these three models clearly suggest that the new method speeds up the computation substantially, and it can be used as a powerful method to estimate the unknown parameters of large-scale network models.

This work shows a common weakness of the inference algorithms, namely a small simulation error does not mean that the generated sample has good accuracy with respect to the exact parameter. This phenomenon can be observed in the third test system model. One of the potential reasons may be that, when the number of model parameters is large, different samples with a variety of values may realize similar simulations with good accuracy.

In this work, we used the tolerance threshold as the key criterion to accept or reject particles. However, for the third example, when the error threshold was small enough, namely the generated solution and observation data almost coincided, it was still difficult to obtain accurate estimated parameters. The possible reason may be that, for a model with a large number of parameters, we may obtain similar simulations by using quite different sets of model parameters. Another issue is how to select the tolerance threshold to optimize and balance the accuracy and efficiency of the inference algorithm. Furthermore, the designed algorithm is an adaptive process to cluster cells and infer the models for cells in the same cluster. This technique may be applied to other inference methods. All of these issues will be interesting topics for further research.