Scalable inference of cell differentiation networks in gene therapy clonal tracking studies of haematopoiesis

Abstract

Motivation: Investigating cell differentiation under a genetic disorder offers the potential for improving current gene therapy strategies. Clonal tracking provides a basis for mathematical modelling of population stem cell dynamics that sustain the blood cell formation, a process known as haematopoiesis. However, many clonal tracking protocols rely on a subset of cell types for the characterization of the stem cell output, and the data generated are subject to measurement errors and noise.

Results: We propose a stochastic framework to infer dynamic models of cell differentiation from clonal tracking data. A state-space formulation combines a stochastic quasi-reaction network, describing cell differentiation, with a Gaussian measurement model accounting for data errors and noise. We developed an inference algorithm based on an extended Kalman filter, a nonlinear optimization, and a Rauch-Tung-Striebel smoother. Simulations show that our proposed method outperforms the state-of-the-art and scales to complex structures of cell differentiations in terms of nodes size and network depth. The application of our method to five in vivo gene therapy studies reveals different dynamics of cell differentiation. Our tool can provide statistical support to biologists and clinicians to better understand cell differentiation and haematopoietic reconstitution after a gene therapy treatment. The equations of the state-space model can be modified to infer other dynamics besides cell differentiation.

Availability and implementation: The stochastic framework is implemented in the R package Karen which is available for download at https://cran.r-project.org/package=Karen. The code that supports the findings of this study is openly available at https://github.com/delcore-luca/CellDifferentiationNetworks.

Figure 1.

Analysis’ flowchart: a clonal tracking dataset and the biochemical reactions (top-right) are the input of our framework Karen (left). It mainly consists of three parts: a filtering step, a maximum likelihood step, and a smoothing step. These steps are iterated until convergence is reached. The inferred cell differentiation network is returned (bottom-right).

Figure 2.

Graphical representation of the cell differentiation networks proposed as candidate models in the in vivo studies (a-d). Grey and white nodes represent unobserved and observed cell types. Light-grey nodes are the cell lineages whose data were collected in all studies except for the genotoxicity one. Arrows represent cell duplication (green), cell death (red), and cell differentiation (blue).

Figure 3.

Graphical representation of the cell differentiation networks used in the in-silico studies (a-d). Grey and white nodes represent unobserved and observed cell types. Arrows represent cell duplication (green), cell death (red), and cell differentiation (blue).

Figure 4.

For each comparative synthetic study with observed (left) and systematically missing (right) progenitors HSC, P1, and P2: boxplots (y-axis) of the estimated parameters divided by the true ones for each reaction rate (x-axis) obtained from each method (colours), across all simulations, under a fraction $\zeta =90\%$ of false negative errors (top), a sampling frequency τ = 4 (middle), and a measurement noise generated by ${\rho }_0={\rho }_1=10$ (bottom).

Figure 5.

(a) For each true data generating process (row) and for each candidate model (column), the simulated process (empty dots), the synthetic data (full dots), the estimated smoothing moments (lines) of a single clone for each cell type (colours), and the median AIC across all simulations used to evaluate model misspecification. (b) Boxplots (y-axis) of the ratio between the estimated and true parameters for each reaction (x-axis) of the cell differentiation network of Fig. 3d used in the scalability study. (c) Estimated smoothing moments (lines), the true Markov states (empty dots), and the synthetic data (full dots) of one clone for each cell type (colours) for a single simulation of the cell differentiation network of Fig. 3d used in the scalability study.

Figure 6.

Inferred cell differentiation networks having the lowest AIC, as defined by Equation (18), for the genotoxicity study for the comparison of the two viral vectors PGK and SFV, the rhesus macaque (RM) study, and the clinical trials WAS, β0βE and βSβS. Each arrow is weighted and coloured according to the transition probabilities estimated with Equation (15).