Mathematical deconvolution of CAR T-cell proliferation and exhaustion from real-time killing assay data

Chimeric antigen receptor (CAR) T-cell therapy has shown promise in the treatment of haematological cancers and is currently being investigated for solid tumours, including high-grade glioma brain tumours. There is a desperate need to quantitatively study the factors that contribute to the efficacy of CAR T-cell therapy in solid tumours. In this work, we use a mathematical model of predator–prey dynamics to explore the kinetics of CAR T-cell killing in glioma: the Chimeric Antigen Receptor T-cell treatment Response in GliOma (CARRGO) model. The model includes rates of cancer cell proliferation, CAR T-cell killing, proliferation, exhaustion, and persistence. We use patient-derived and engineered cancer cell lines with an in vitro real-time cell analyser to parametrize the CARRGO model. We observe that CAR T-cell dose correlates inversely with the killing rate and correlates directly with the net rate of proliferation and exhaustion. This suggests that at a lower dose of CAR T-cells, individual T-cells kill more cancer cells but become more exhausted when compared with higher doses. Furthermore, the exhaustion rate was observed to increase significantly with tumour growth rate and was dependent on level of antigen expression. The CARRGO model highlights nonlinear dynamics involved in CAR T-cell therapy and provides novel insights into the kinetics of CAR T-cell killing. The model suggests that CAR T-cell treatment may be tailored to individual tumour characteristics including tumour growth rate and antigen level to maximize therapeutic benefit.

. CAR T-cell therapy is a specific type of immunotherapy where T-cells are genetically modified to recognize a tumour antigen thereby specifically redirecting T-cell cytolytic activity. Inspired by the success of CAR T-cell therapy in liquid tumours, there has been great interest in expanding the use of CAR T-cells for the treatment of solid tumours, such as glioblastoma (GBM), a highly aggressive form of primary brain cancer. Several clinical trials using CAR T-cells to treat GBM have been initiated all over the world [3][4][5][6]. At this early stage of clinical development, CAR T-cells offer much promise in solid tumours. However, the diversity of current clinical trials employing varying types of CARs for different solid tumours, target patient populations, preconditioning regimes and cell origins (autologous versus allogeneic) presents a significant challenge in identifying which aspects of a given CAR T-cell treatment protocol are most critical for its effectiveness. An additional critical challenge for CAR Tcell therapy is the potential for transient progression, where the cancer appears to progress before eventually responding to the treatment [7,8].
In order to address these challenges in CAR T-cell therapy for solid tumours, we endeavoured to study the kinetics of CAR T-cell killing with an in vitro system and a mathematical model. Mathematical models are useful to describe, quantify and predict multifaceted behaviour of complex systems, such as interactions between cells. A mathematical model is a formalized method to hypothesize systems dynamics, and yield solutions that predict the system's behaviour with a given set of parameters and initial conditions. Mathematical models can be versatile and tested with clinical data which may be obtained in vivo from non-invasive imaging [9][10][11] and the models can be refined when additional information about the system becomes available. Many mathematical models have been developed to understand tumour progression to guide refinement of cancer therapy regimens [12][13][14].
As CAR T-cell therapy is a newly advanced treatment modality, relatively few studies have used computational modelling to understand and improve this cell-based therapy. Recently, computational models have been developed to investigate cytokine release syndrome for toxicity management [15][16][17], effect of cytokine release syndrome on CAR T-cell proliferation [18], mechanisms of CAR T-cell activation [19,20], and dosing strategies [21]. However, it remains an open challenge how to use mathematical modelling to study and ultimately predict dynamics of CAR Tcell mediated cancer cell killing with respect to CAR T-cell dose, donor-dependent T-cell differences, cancer cell proliferation, target antigen expression, and how these factors contribute to the overall effectiveness of CAR T-cell therapy.
Based upon our pre-clinical and clinical experience with our well-characterized IL13Rα2-targeted CAR T-cell therapy for recurrent GBM [22,23], we have identified several factors which contribute to the effectiveness of CAR T-cells, namely: rates of proliferation, exhaustion, persistence and target cell killing. To study these various facets of CAR T-cell killing kinetics, we modelled the dynamics between cancer cells and CAR T-cells as a predator-prey system with a mathematical model we call CARRGO: Chimeric Antigen Receptor T-cell treatment Response in GliOma. We use a real-time cell analyser experimental system to estimate parameters of the mathematical model and then apply the model to in vivo human data. The long-term aim of this work is to develop a model which could be used to predict and eventually to optimize response to CAR T-cell therapy.

Methods
The CARRGO mathematical model is a variation on the classic Lotka-Volterra [24,25]  where X represents the density of cancer cells, Y is the density of CAR T-cells, ρ is the net growth rate of cancer cells, K is the cancer cell carrying capacity, k 1 is the killing rate of the CAR T-cells, k 2 is the net rate of proliferation including exhaustion of CAR T-cells when encountered by a cancer cell and u is the death rate of CAR T-cells. The parameters ρ, K, k 1 , u are constants and assumed to be non-negative except for k 2 which can be either positive or negative (table 1). A positive value of k 2 indicates an increased rate of CAR T-cell proliferation when stimulated by interaction with a cancer cell. A negative value of k 2 indicates exhaustion or limited activation of CAR T-cells resulting from interaction with a cancer cell. Exhaustion and hypoactivation of CAR T-cells are combined into a single value and are not modelled individually. We chose to model the net number of cancer cells and simple interactions between cancer cells and CAR T-cells because the output data from the culture system are limited to cell number over time. We therefore are only able to infer dynamics at this scale and dimension (i.e. number of cells and time). Moreover, we performed a system identifiability analysis to demonstrate the parameters of the model can be uniquely determined from the data in this experiment (see electronic supplementary material, methods) [26][27][28][29]. Future studies may examine more complex dynamics such as individual cell antigen levels, heterogeneity, resistant and sensitive sub-populations, repeated treatments, etc. with other experimental designs which directly measure these features.

Model assumptions
The CARRGO model treats cancer cell-CAR T-cell dynamics in this experimental condition as a closed predator-prey system. The model assumes (1) the populations are well mixed, (2) cancer cell growth is limited by space and nutrients (culture media) in the in vitro culture system and therefore grow logistically, (3) CAR T-cells kill cancer cells when they interact via the law of mass action, (4) the CAR T-cell killing rate does not explicitly assume a dependence on antigen density, (5) CAR T-cells may be stimulated to proliferate or to undergo loss of effector function-defined as exhaustion-upon contact with a cancer cell [30], and (6) the CAR T-cell death rate is independent of cancer cell density. We chose the logistic growth model for the cancer cell population because the fixed growth rate and carrying capacity parameters were the biological quantities of interest when comparing CAR T-cell killing kinetics across cell lines. Witzel et al. compared several sigmoidal growth laws including logistic, Gompertz and Richards, and showed that all these models can be fitted equally well to this form of experimental data [31]. Data supporting our model assumptions are given in electronic supplementary material, figures S1 and S2.

Dynamical system analysis of the CARRGO model
Closed form solutions cannot be obtained for the relatively simple CARRGO model. To study the possible dynamics of the CARRGO model, we perform classical dynamical system analysis. Detailed mathematical analysis of this model can be found in several textbooks on dynamical systems [25,32]. In the interest of informing the reader, we briefly summarize the main points here. We begin by (1) scaling (non-dimensionalizing) the variables in the system and then (2) identify stationary points and classify their stability and finally (3) interpret the stationary points and system dynamics in terms of the initial numbers of cancer cells and CAR T-cells. First, we scale the variables in the CARRGO model to obtain a model without physical units in order to study the intrinsic dynamics of the system. We scale time, the cancer cell and CAR T-cell populations as t ¼ tr, y ¼ ðk 1 =rÞY, x ¼ X=K: These variables are substituted into the CARRGO model (equations (3.1) and (3.2)) to obtain the scaled dimensionless system The steady-state solutions of this system (equations   cancer cells. In this case, the fixed points (0,y) and (K,0) are the same as in case 1. However, the point (K,0) is now a stable attractor state, corresponding to the extinction of CAR T-cells and the cancer cells eventually growing to carrying capacity. Again, cancer cell elimination (0,y) is stable only if the CAR T-cell population is larger than the ratio of cancer cell proliferation and CAR T-cell killing rate ðY . r=k 1 Þ. If the initial CAR T-cell population satisfies the inequality Y . ðr=k 1 Þ À ðrX=Kk 1 Þ, the CARRGO model predicts a transient regression of cancer cells before eventual rapid progression (grey region). This case is a failure of CAR T-cell treatment (figure 1b). This case results in an oscillating behaviour of an increase in cancer cells corresponding to tumour progression followed by a decrease in tumour cells corresponding to treatment response. The transient and oscillatory nature of these dynamics may be interpreted as a 'pseudo'-failure and 'pseudo'-response to the therapy. We note that cancer progression and treatment occur on finite and sometimes small time scales and therefore oscillatory dynamics may not be observed in vivo due to insufficient time to observe these changes.

Cell lines
Low-passage primary brain tumour (PBT) lines were derived from GBM patients that had undergone tumour resections at City of Hope as previously described [23,33]. Fibrosarcoma line HT1080 was obtained from the American Tissue Culture Collection (ATCC) and maintained according to recommendations. The cell line PBT030 endogenously expresses high level of IL13Rα2. HT1080 and PBT138 do not express IL13Rα2 and were lentivirally engineered to express varied levels based on different promoter strengths to investigate the relationship between killing kinetics and antigen expression level: high (greater than 70% + ) driven by the EF1α promoter, medium (between 40% + and 70% + ) driven by the PGK promoter, low (less than 20% + ) driven by the attenuated PGK100 promoter [34,35]. These cell lines are denoted with H, M, L, respectively (e.g. HT1080-H). These tumour cell lines were selected because they differ in aggressiveness (proliferation rates) and antigen expression levels (endogenous or engineered).
CAR T-cells were derived from healthy donor CD62 L+ CD45RO+ central memory T-cell population and lentivirally transduced with second-generation IL13Rα2-targeting CARs: IL13BBζ or IL1328ζ [23,33,36,37]. Transduced product was enriched for CAR and expanded in X-Vivo media with 10% fetal bovine serum until 17 days in culture and cryopreserved. Non-transduced T-cells expanded under the same condition were used as mock control.

Experimental design
Real-time monitoring of cancer cell growth was performed by using xCELLigence cell analyser system [38]. This system uses electrical impedance to non-invasively quantify adherent cell density with a dimensionless number referred to as cell-index (CI). The cell-index read-out is strongly positively correlated with the number of cells in the well (r 2 > 0.9) and can be used as a measure of cell number [31]. We therefore report CARRGO parameter values in units CI which can be translated into units per cell based on the linear relation. Real-time cytotoxicity assay was performed using xCELLigence system in disposable 96 well E-Plates. Prior to seeding, tumour cells were enzymatically single-celled and seeded at 25 × 10 3 , 12.5 × 10 3 or 2 × 10 3 cells per well depending on the cell line. Cells were either left untreated (triplicates per cell line) or treated with CAR T-cells at effector to target ratios (E : T) of 1 : 5, 1 : 10 and 1 : 20. CAR T-cells were added to the wells about 24 h after cancer cell seeding. Growth curves were recorded over 4 days with temporal resolution of 15 min (figure 2). Each cell line was treated with three IL13Rα2-targeted CAR T-cells: BBζ, 28ζ and mock. At the end of the experiment, flow cytometry was performed to measure the residual CAR T-cells, cancer cells and IL13Rα2 expression level. The details of cancer cell seeding and effector to target ratios used for the experiments are given in electronic supplementary material, table S1. The cancer cell dynamics of all the wells of 96 well E-plate for all cell lines are given in electronic supplementary material, figure S3.

CARRGO model fitting to experimental data
The first 24 h of the time series describes the process of cell attachment to the bottom of the plate (figure 2). The spatial process of cell adhesion and spreading in the well can be modelled as a reaction-diffusion process, described in electronic supplementary material, figure S4. Since we are interested in cell growth kinetics, we omitted the data from the first 24 h during the attachment process. The final time point is determined when the cells reach confluency which varies by cell line. Observed time of confluency for the PBT cell line was around 120 h from the time of seeding while HT1080 reached confluency within 80 h. The data points from CAR T-cell administration (24 h after seeding) up to 80% of maximum CI value (confluency) were used for model fitting to estimate the parameters. At greater than 80% of the maximum CI, the linear relationship between CI and cell number no longer holds [31].
Cancer cell net growth rate r and carrying capacity K (equation (3.1)) were computed by fitting logistic growth to untreated cancer cell time-series data. The CAR T-cell killing rate k 1 , growth rate k 2 and death rate u were computed by fitting the solution of the CARRGO model (equations (3.1) and (3.2)) to treated cancer cell time-series data by minimizing the root mean square error. Linear regression was used to determine the quality of model fitting (R 2 ). All optimization computations were performed in Matlab with fmincon.

CARRGO model fitting to human data
A patient with recurrent glioma received CAR T-cells engineered for IL13Rα2 and showed complete tumour regression, which was published as a brief report by Brown et al. in 2016 [23]. We retrospectively collected the magnetic resonance imaging (MRI) data of this patient and calculated tumour volumes to be used to fit the CARRGO model. Three lesions were selected using the lesion labelling reported in Brown et al. [23]: lesions T6, T7 which responded to IL13Rα2 targeted CAR T-cells and lesion T9 which was a lesion that appeared later on which did not respond to the therapy. Tumour volumes for each lesion were estimated by manual segmentation of contrast-enhancing lesions from T1-weighted post-contrast MRIs. The number of cancer cells (CC) was estimated by calculating CC = [tumour volume (µm 3 )]/ [GBM cell size (µm 3 )] with the average cell diameter assumed to be 20 µm [39]. The volume of a spherical cell is then given by V ¼ ð4p=3Þð10 mmÞ 3 . These relationships were used to estimate the total number of cancer cells in a tumour volume. The tumour growth rate (ρ) for lesions T6 and T7 was computed from two subsequent imaging time points following the first appearance of the lesion on MRI. Lesion T9 underwent surgical resection prior to CAR T-cell administration. The tumour growth rate for lesion T9 was computed from the pre-surgical MRI data. CAR T-cells were administered first directly into the tumour tissue and subsequently into the cerebrospinal fluid via royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20190734 the intraventricular injection. Because the CAR T-cells migrated to several tumour foci in the patient, we assumed a small fraction (5-10%) of the infused dose reached each individual lesion at each infusion. The CARRGO model was fitted to the timeseries MRI-derived tumour volume data by minimizing the root mean square error for MRIs before and during CAR T-cell treatment to compute the rates of CAR T-cell killing k 1 , exhaustion k 2 and death θ. For lesion T9, the CARRGO model was fitted to the MRI data following CAR T-cell administration which occurred after the partial surgical resection.

Model/data fitting to in vitro data
A high goodness of fit of the CARRGO model to the xCELLigence data was observed across all cell lines (R 2 ¼ 0:93 + 0:1, figure 3; electronic supplementary material, figure S3). To investigate the sensitivity of our model fitting to sampling frequency, we down-sampled the data by taking time intervals of 2, 5 and 10 h. No significant variation was observed in the model parameters k 1 , k 2 and u to the down-sampled data (repeated measure ANOVA p > 0.1) (electronic supplementary material, figure S5 and S6). We consistently observed very small values of the CAR T-cell death rate ðu , 10 À3 Þ. Uniqueness of the parameters was tested by choosing 100 different combinations of values of the parameters across several orders of magnitude for the model fitting optimization procedure. We found that if the optimization converged, it converged to unique values of the parameters, which is a direct consequence of the identifiability analysis of the model and minimum number of points required to resolve the model (see electronic supplementary material, S1, figure S7 and movie S1).

CAR T-cell dose-dependent dynamics
We examined the effect of varying the effector to target ratio, i.e. CAR T-cell dose, for all cell lines. The CAR T-cell death rate parameter was found to be very small (u , 10 À3 day −1 ) for all cancer cell lines and all CAR T-cells and doses. The killing rate parameter k 1 shows a negative correlation while k 2 shows positive correlation with respect to the CAR T-cell dose. The parameter k 2 was negative only for tumour line HT1080-H which indicates the exhaustion rate being much stronger than the proliferation rate of the CAR T-cells.

Relating k 1 , k 2 with tumour growth rate and antigen expression
Tumour growth rate r varies significantly ( p < 0.01) among different cell lines and with antigen expression level (see electronic supplementary material S2 and figure S11a). To investigate the relationship between tumour growth rate and CAR T-cell killing k 1 and exhaustion k 2 , we evaluated cell lines with antigen levels greater than 80% and treated with BBζ CAR T-cells at an effector to target ratio of 1 : 5. No significant correlation was found between the cancer cell proliferation rate r and killing rate (k 1 ) (electronic supplementary material, figure S11b). However, the exhaustion rate k 2 is significantly correlated with tumour growth rate (electronic supplementary material, figure S11c) with Pearson correlation coefficient r ¼ À0:9, p < 0.001. Similar results were observed for the cells treated with 28ζ IL13Rα2-CARs. Figure 5 shows the density of IL13Rα2 level on cancer cell surface and its relation to CAR T-cell killing for cell line HT1080-H and PBT138-H. We observed that k 1 shows a decreasing trend from medium to high antigen level (figure 5c) suggesting that high levels of antigen expression may not result in faster rates of CAR T-cell killing. The rate constant k 2 increases from low to medium antigen expression and plateaus with high levels ( figure 5d ). This suggests limited activation of CAR T-cell at lower antigen expression and exhaustion rate from medium to high antigen may not change significantly and may be the result of over-activation of the CAR T-cells.

CARRGO model applied to in vivo human data
To translate the in vitro dynamics of the model to clinical data [23], we fit the CARRGO model to MRI-derived tumour volumes during CAR T-cell treatment ( figure 6). The CARRGO model is able to fit the tumour growth dynamics quite accurately for lesions T6 and T7 with the same set of parameters k 1 ¼ 6 Â 10 À9 (day −1 cell −1 ), k 2 ¼ 0:3 Â 10 À10 (day −1 cell −1 ), u ¼ 0:1 Â 10 À5 (day −1 ) and lesion T9 with k 1 ¼ 9 Â 10 À8 (day −1 cell −1 ), k 2 ¼ À2 Â 10 À13 (day −1 cell −1 ), u ¼ 5 Â 10 À5 (day −1 ). In the case of lesion T9, although the CARRGO model is consistent with the overall tumour dynamics, it does not fit the later time points following CAR T-cell treatment well. This is because lesion T9 received radiation treatment between day 200 and 300, which is not included in the CARRGO model. We note the negative correlation between the tumour growth rate (r ¼ 0:06, 0:07 and 0.2 day −1 for T6, T7 and T9, respectively) with the CAR T-cell exhaustion rate k 2 in the patient data, which is consistent with that observed in the experimental data (electronic supplementary material, figure S11c). We remark that the parameters k 1 and k 2 are of the order of O(10 À13 ), which appear to be very small; however, these parameters are scaled by the carrying capacity K in units of cells, which is of order Oð10 9 Þ. Therefore, these parameter values are comparable with the in vitro data when scaled relative to the carrying capacity (figure 4).

Discussion
The CARRGO model is a simple representation of cancer cell-T-cell interactions. We developed the CARRGO model with the aim of understanding CAR T-cell efficacy in terms of rates of killing, proliferation, exhaustion, and persistence with a real-time cell analyser in a simple, controlled in vitro system. The CARRGO model fitted remarkably well to the highly temporally resolved experimental data and as well as to data derived from a patient treated with IL13Rα2 BBζ CAR T-cells. Although the predator-prey mathematical model formalism has been widely used in a number of biological settings, the novelty of this model is in the application to a novel form of cancer therapy with a high temporal resolution cell monitoring experimental design which provides nearly continuous data on killing kinetics.
With the CARRGO model we show that the rate of cancer cell killing by CAR T-cells is inversely related to the CAR T-cell dose. With a fixed number of cancer cells as an initial condition, as the number of CAR T-cells increases (dose), any individual T-cell will encounter fewer number of cancer cells to kill, indicating that increasing dose does not result in a maximal rate of killing on a per T-cell basis. For example, the PBT138 cell line shows complete killing within 80 and 100 h for effector to target ratios 1 : 5 and 1 : 10, respectively. This result suggests that a lower dose of CAR T-cells may change the time to complete cancer cell killing but shows the same overall cancer cell killing effectiveness. Moreover, we observed that k 2 positively correlated with CAR T-cell dose. Because the parameter k 2 is a net measure of CAR T-cell proliferation and exhaustion or lack of activation and the T-cell proliferation is not dose-dependent [18], the trend observed in k 2 with dose is dominated by the exhaustion rate: the higher the dose, the lower the exhaustion rate, resulting in an increased value of k 2 . The death rate of CAR T-cells was very small when compared with the cancer cell proliferation rate for all conditions. This is likely due to the short  the clinical phenomenon of pseudo-progression, in which the cancer is seen to progress during therapy before eventually responding [7,8]. Identifying characteristics of the patient and the CAR T-cells which may result in pseudo-progression could have a profound effect on interpretation of these dynamics observed in the clinic. Interestingly, we found k 1 decreases and k 2 plateaued from medium antigen level to higher level of antigen expression. One of the possible explanations of this behaviour could be the antigen density is more heterogeneous in the higher antigen level cell population when compared with medium and low antigen levels (figure 5a). More heterogeneity in the density of antigen expression intensity in the cancer cells within the initial population may cause clustering of CAR T-cells resulting in their exhaustion [20,40]. Another confounding factor can be the dependence of the detected antigen signal intensity on both the number of antigen-positive tumour cells and their individual antigen expression intensity. Elucidating these factors individually can better tune the model parameters and the prediction of the tumour response dynamics. However, more studies are required to examine effect of cell and population level antigen density on CAR T-cell killing kinetics.
There are some important limitations to consider with this model and experimental system. Perhaps the most obvious is that the in vitro system is not a model for the human immune system or tumour microenvironment. It does not include cytokines, stromal cells or additional immune cells such as myeloid cells which contribute to CAR T-cell activity in vivo. Another limitation is the assumption that the populations are well mixed. In practice, this assumption may depend on the route of CAR T-cell administration, with intracavitary and intraventricular injections potentially resulting in spatially heterogeneous densities of CAR T-cells or cancer cells, although methods to assess the distribution of CAR T-cells in vivo remain an open challenge [41,42]. To address this limitation, the well-mixed assumption may be relaxed and CAR T-cell killing dynamics interrogated with spatial or agent-based models [43]. Another limitation is with regard to the experimental system: the change in electrical impedance measured by the cell index does not differentiate cell detachment from cell killing. This is only a minor consideration as the cell lines used are very adherent to the plate and were not observed to detach. Finally, the experimental system does not directly measure dynamics of CAR T-cells. However, our model is initialized with known numbers of cancer cells and CAR T-cells and the model-predicted cancer cell to CAR T-cell ratio at the experimental endpoint was validated with flow cytometry, giving confidence to our model predictions and parameter estimates. To address this limitation, the CARRGO model CAR T-cell dynamics can be validated by labelling the CAR T-cells and directly measuring their dynamics with live cell imagingbased methods [44]. Despite these limitations, the CARRGO model succeeded in revealing nonlinear dynamics, quantifying kinetics of killing, and generating hypotheses which may be tested in other in vitro systems, and other computational or in vivo models.
An interesting application of this model would be to adapt the system to evaluate differences between donor T-cells-either autologous or allogeneic. The intrinsic fitness of the T-cell used for CAR T-cell manufacturing is known to be a critical differentiator between responding and non-responding patients [45], and this model may be able to predict T-cell products with high proliferative potential versus products more prone to exhaustion. While our studies used an allogeneic co-culture platform for the predator-prey modelling of effector activity, we expect these findings to be translatable to autologous CAR T-cell therapies, as CAR T-cells recognize targets in an MHC-independent manner. Future studies could specifically evaluate allogeneic and autologous donor-dependent differences to assess the fit of the mathematical model, and guide potential adjustments to the model. The ability to differentiate key quality attributes of the therapeutic product could potentially provide powerful predictions of CAR T-cell product efficacy and facilitate patient treatment management.
In summary, CAR T-cells have shown promise in hematologic malignancies and are being actively investigated in solid tumours. We aimed to use mathematical modelling to investigate factors which contribute to the kinetics of CAR T-cell mediated cancer cell killing in a simple isolated in vitro system. We were able to fit the CARRGO model to in vitro and in vivo human data with remarkable accuracy. We demonstrated that we can consistently and reproducibly estimate rate constants in the CARRGO model and investigate their dependence on CAR T-cell dose and antigen expression levels. The CARRGO model may be combined with other mathematical models which estimate cancer cell growth and proliferation rates non-invasively with MRI data [9,11,46] to produce a fine-tuned and benchmarked suite of mathematical models, which may aide in optimization of dosing and scheduling of CAR T-cells for greater individualized and personalized therapy.