Data driven model discovery and interpretation for CAR T-cell killing using sparse identification and latent variables

Abstract

In the development of cell-based cancer therapies, quantitative mathematical models of cellular interactions are instrumental in understanding treatment efficacy. Efforts to validate and interpret mathematical models of cancer cell growth and death hinge first on proposing a precise mathematical model, then analyzing experimental data in the context of the chosen model. In this work, we present the first application of the sparse identification of non-linear dynamics (SINDy) algorithm to a real biological system in order discover cell-cell interaction dynamics in in vitro experimental data, using chimeric antigen receptor (CAR) T-cells and patient-derived glioblastoma cells. By combining the techniques of latent variable analysis and SINDy, we infer key aspects of the interaction dynamics of CAR T-cell populations and cancer. Importantly, we show how the model terms can be interpreted biologically in relation to different CAR T-cell functional responses, single or double CAR T-cell-cancer cell binding models, and density-dependent growth dynamics in either of the CAR T-cell or cancer cell populations. We show how this data-driven model-discovery based approach provides unique insight into CAR T-cell dynamics when compared to an established model-first approach. These results demonstrate the potential for SINDy to improve the implementation and efficacy of CAR T-cell therapy in the clinic through an improved understanding of CAR T-cell dynamics.

Keywords: CAR T-cells; SINDy; allee effect; antigen binding; cell therapy; dynamical systems; glioblastoma; latent variables.

Figure 1

Diagram of experimental procedure highlighting use of microelectrode plates in an xCELLigence cell analyzer system and sample Cell Index (CI) measurements for control and treatment groups (E:T = 1:4). This system utilizes real-time voltage measurements to determine CI values representative of the adherent cancer cell population as a function of time. CAR T-cells are added following 24 hours of cancer cell expansion and attachment. After 6-8 days of monitoring the cancer cell growth and death dynamics, cells are harvested and enumerated using flow cytometry.

Figure 2

Conceptual graphs of population size (in cell index - CI) versus time (in hours - hrs) for the three growth models presented in Eqs. (3)-(5): logistic growth (A) weak Allee effect (B) strong Allee effect (C). Model parameter values are: $\rho =0.75$ hrs^-1, $K=10$ CI, $A=5$ CI, and $B=5$ CI. Colors correspond to different initial cell populations, which are the same for each model presented (blue = 12 CI, orange = 8 CI, green = 4 CI, red = 1 CI).

Figure 3

(A) Compartmental model for single and double CAR T-cell-cancer cell binding. Expressions for how rate constants ( $k_i^{(j)}$ ) contribute to the growth or death of the cancer cell and CAR T-cell populations are presented in Eqs. (9)-(14). See (30) for further development and analysis of the cell binding model. (B) Graphs of binding rate versus CAR T-cell population for the single binding, double binding, and effective double binding models in Eqs. (9)-(12), (16), and (18). Model parameters for antigen bindings are: $a=20$ CI^-2 hrs^-2 and $h=16$ CI^-1 hrs^-1 for single binding; $a=20$ CI^-2 hrs^-2, $b=5$ CI^-3 hrs^-2, $h=16$ CI^-1 hrs^-1, and $k=2$ CI^-2 hrs^-1 for double binding; and $a=20$ CI^-2 hrs^-2, $b=2.75$ CI^-3 hrs^-2, $h=16$ CI^-1 hrs^-1, and $k=2$ CI^-2 hrs^-1 for effective double binding. These parameter values were chosen to highlight how well the effective double binding model can approximate both the single and double binding models at low CAR T-cell population values, $y<1$ CI. Note that since the original double binding model in this scenario is concave-up, the effective double binding model parameters should be chosen to match concavity. This requirement sets a positivity constraint on the quadratic term in Eqs. (16) and (18). (C) Graphs of CAR T-cell response rates versus cancer cell population for different functional response models. Model parameters for functional responses are: $p=6/5$ CI^-1 hrs^-1 for Type I; $p=$ CI^-1 and $g=5$ CI for Types II and III. Note overlap of Types I and II functional responses for $x<1$ CI, and distinct differences in concavity between Types II (negative) and III (positive) for $x<2$ CI. These characteristics correspond to Type I and Type II functional responses being indistinguishable at low cancer cell populations, and Type II and Type III being differentiated by fast-then-slow response rates (Type II) versus slow-then-fast response rates (Type III).

Figure 4

(A-C) Latent variable analysis results for first of two experimental replicates each E:T ratio examined. Presented are the cancer cell index measurements from the xCELLigence machine in red, overlaid with the splined measurements for the cancer cells in black; the two endpoint measurements for the CAR T-cell levels enumerated by flow cytometry in black, with the CAR T-cell population trajectory as determined by latent variable analysis in yellow, overlaid with the splined CAR T-cell trajectory in blue. Note that despite the CAR T-cell populations being measured with flow cytometry, we have converted levels to units of Cell Index for ease of comparison with the cancer cells, using a conversion factor of 1 CI ≈ 10,000 cells. (D-F) Predicted trajectories of discovered models compared to splined measurements of cancer cells and CAR T-cells for same data presented in (A-C). Splined cancer cell and CAR T-cell measurements are in black and blue, respectively. Predicted trajectories for cancer cells are the red dot-dashed lines, while the CAR T-cells are the purple dot-dashed lines. To examine stability of SINDy-discovered models, both simulations and forward predictions are presented to show steady-state behavior. Note that the best fits between predictions and measurements occur in the high E:T scenario, where assumptions made regarding treatment success and low cancer cell populations in determining model candidate terms are best adhered. As the E:T ratios get smaller, increasing deviation between discovered model predictions and splined measurements can be qualitatively observed. This is likely due to weakening of assumptions of treatment success and low cancer cell populations associated with the low E:T conditions. See Supplemental Material Figure S2 for equivalent latent variable analysis results and SINDy-predicted trajectories for the second set of experimental replicates.

Figure 5

Predictions of cell trajectories for E:T ratios of 1:4, 1:8, and 1:20 from CARRGO model (blue) and SINDy model (red). Model fits for both CARRGO and SINDy were performed using Levenberg-Marquadt Optimization (LMO) on data aggregated across experimental replicates. Initial LMO parameter value guesses were determined by parameter values from SINDy or from published CARRGO model values. Data points represent the mean of all experimental replicates, while error bars represent the ranges across replicates. Of note are the differences in CARRGO and SINDy model predictions for the final CAR T-cell values compared to measurements, and the notable difference in when the maximum CAR T-cell population is reached between CARRGO and SINDy models. Note that experimental measurements have been down-sampled to 25% to allow for visualization.