Predicting Cytotoxic T cell Age from Multivariate Analysis of Static and Dynamic Biomarkers

Adoptive T-cell transfer therapy relies upon in vitro expansion of autologous cytotoxic T cells that are capable of tumor recognition. The success of this cell-based therapy depends on the specificity and responsiveness of the T cell clones before transfer. During ex vivo expansion, CD8+ T cells present signs of replicative senescence and loss of function. The transfer of nonresponsive senescent T cells is a major bottleneck for the success of adoptive T-cell transfer therapy. Quantitative methods for assessing cellular age and responsiveness will facilitate the development of appropriate cell expansion and selection protocols. Although several biomarkers of lymphocyte senescence have been identified, these proteins in isolation are not sufficient to determine the age-dependent responsiveness of T cells. We have developed a multivariate model capable of extracting combinations of markers that are the most informative to predict cellular age. To acquire signaling information with high temporal resolution, we designed a microfluidic chip enabling parallel lysis and fixation of stimulated cell samples on-chip. The acquisition of 25 static biomarkers and 48 dynamic signaling measurements at different days in culture, integrating single-cell and population based information, allowed the multivariate regression model to accurately predict CD8+ T-cell age. From surface marker expression and early phosphorylation events following T-cell receptor stimulation, the model successfully predicts days in culture and number of population doublings with R2=0.91 and 0.98, respectively. Furthermore, we found that impairment of early signaling events following T cell receptor stimulation because of long term culture allows prediction of costimulatory molecules CD28 and CD27 expression levels and the number of population divisions in culture from a limited subset of signaling proteins. The multivariate analysis highlights the information content of both averaged biomarker values and heterogeneity metrics for prediction of cellular age within a T cell population.


Summary
Adoptive T cell transfer therapy relies upon in vitro expansion of autologous cytotoxic T cells that are capable of tumor recognition. The success of this cell-based therapy depends on the specificity and responsiveness of the T cell clones before transfer. During ex-vivo expansion, CD8+ T cells present signs of replicative senescence and loss of function. The transfer of nonresponsive senescent T cells is a major bottleneck for the success of adoptive T cell transfer therapy. Quantitative methods for assessing cellular age and responsiveness will facilitate the development of appropriate cell expansion and selection protocols. Although several biomarkers of lymphocyte senescence have been identified, these proteins in isolation are not sufficient to determine the age-dependent responsiveness of T cells. We have developed a multivariate model capable of extracting combinations of markers that are the most informative to predict cellular age. To acquire signaling information with high temporal resolution, we designed a microfluidic chip enabling parallel lysis and fixation of stimulated cell samples on-chip. The acquisition of 25 static biomarkers and 48 dynamic signaling measurements at different days in culture, integrating single-cell and population based information, allowed the multivariate regression model to accurately predict CD8+ T cell age.
From surface marker expression and early phosphorylation events following T cell receptor stimulation, the model successfully predicts days in culture and number of population doublings with R 2 = 0.91 and 0.98, respectively. Furthermore, we found that impairment of early signaling events following T cell receptor stimulation due to long term culture allows prediction of co-stimulatory molecules CD28 and CD27 expression levels and the number of population divisions in culture from a limited subset of signaling proteins. The multivariate

Introduction
Immune cell-based therapies hold promise in cancer therapy by harnessing the body's natural defense mechanisms against tumors, while leaving healthy cells unharmed (1,2). Among those therapies, adoptive transfer of T cells has resulted in encouraging clinical trials for treating metastatic melanoma as well as non-Hodgkin's lymphoma, chronic lymphocytic leukemia and neuroblastoma (3)(4)(5). Although cancer cells are less immunogenic than pathogens, the adaptive immune system is able to recognize and eliminate tumor cells. Adoptive therapy with cytotoxic CD8+ T cells (CTLs) relies on the isolation of functional and tumor specific T cells and large in vitro clonal expansion. Once transferred back in the cancer patient, CTLs need to retain tumor specificity and proliferate further in vivo to establish an effective in vivo response and tumor shrinkage. In vivo persistence is a critical factor for elimination of residual or recurring malignant cells. The encouraging results of adoptive transfer therapy could be improved by enhancing the quality of transferred T cells. Cells derived from aged cancer patients have a skewed immune repertoire towards cells that underwent extensive clonal expansion against persistent antigens, resulting in few tumor-specific CTLs (6)(7)(8). Once isolated the tumor infiltrating cells go through a prolonged ex vivo culture process. T cells, as other somatic cells, have a finite clonal lifespan. Extensive ex vivo proliferation and clonal expansion result in T cell differentiation and ultimately replicative senescence (9). To obtain sufficient number of cells before transfer, tumor specific CTLs are activated and undergo several rounds of divisions, resulting in the progressive shortening of telomeres. Chronic antigenic stress and critically short telomere length lead CTLs to enter a state of senescence characterized by functional changes.
Although extensively cultured CTLs retain antigen specificity for the tumor (10), they present Predicting Cytotoxic T cell Age 8 in RPMI 1640 medium with L-glutamine (Sigma-Aldrich) with 10 mM HEPES, 1 mM sodium pyruvate, and 1X MEM nonessential amino acids, and 100 units.mL -1 penicillin/streptomycin (Cellgro) and 10% certified heat-inactivated fetal bovine serum (Sigma-Aldrich). The culture medium was supplemented with 50 U/mL of recombinant IL-2 (Sigma-Aldrich) and Dynabeads® Human T-Activator CD3/CD28 (Invitrogen) at 1:1 bead to cell ratio (kept constant for the entire culture period) for rapid cell expansion (23,24). Cell cultures were checked daily and resuspended in fresh medium when needed. The number of population doublings (PDs) was calculated from the average cell count using the following equation: (1) where n post represents the number of cells counted after expansion and n init represents the number of cells initially seeded.

Microfluidic device fabrication
The two-module device was fabricated using standard soft lithographic techniques (20,25).
Briefly, the modules were molded in poly (dimethylsiloxane) (PDMS) (Dow Corning Sylgard 184, Essex-Brownwell Inc.) from a two-layer SU-8 (Microchem Corp.) master. One layer of 70 mthick SU-8 2050 was spun onto 100 mm silicon wafer, prebaked, and exposed under UV light to define a negative image of the channel system in the resist, following the manufacturer's instructions. After postbaking to crosslink the exposed resist, another layer of 40 m thick SU-8 2020 was spun on top. This layer formed the staggered herringbone arrays (20). After the same prebake and expose process, the wafers were developed using propylene glycol monomethyl

Device operation
A syringe pump (Chemyx Fusion 200 series) controlled the flow to the four inlets at 44 L·min -1 .
Before running an experiment, a solution of 2% BSA in PBS was flown through the device with plugged outlets to pressurize the device and remove any air bubbles. Then the outlets were opened to atmospheric pressure and PBS with 2% BSA was flowed through the device for an additional 15 min. Seven hours before the device operation, 20·10 6 cells were resuspended in fresh medium without IL-2 and anti-CD3/CD28 beads. During device operation, cells were delivered in 1.5 mL of a PBS + 7% w/v dextran solution to match the cell density and avoid cell settling (20). To further avoid cell loss in the syringe, a small cubic magnet with 2-mm edges was inserted with the cell suspension, and intermittently agitated during the experiment. The stimulus consisted of PBS supplemented with 2 g/mL of anti-CD3, clone OKT3 (eBioscience), and 2 g/mL of anti-CD28 (BD Bioscience). Ice-cold freshly-prepared lysis buffer as previously described (20) and a 10% formalin solution (Sigma-Aldrich) were delivered at the inlets of the second module. Cells, stimulus, lysis buffer and fixing solution were flowed for 9 min to allow steady state for the different time points to be reached before the 8 fixed samples and 8 lysates were collected in ice-cold 96-well plates covered with paraffin for 20 additional min.

Signaling measurements
Total protein concentration of the lysates was determined with a BCA assay kit (Pierce

Flow Cytometry analysis
Surface marker protein expression was determined by direct immunostaining. Stimulated Cell cycle analysis was performed on the FlowJo cell cycle analysis platform using the Dean Jett Fox model.

Partial Least Square Modeling
Statistical modeling was performed using the SIMCA-P software (Umetrics). All signals were mean centered and unit variance scaled prior to analysis to allow all variables to be considered equally scaled in principal components (27). The data set was divided into two matrices: Y ε R 18*2 consisting of measures of age in culture (dependent variable block), and X ε R 18*140 , denoting the measured protein phosphorylation signals as well as surface markers, cell morphology and intracellular proteins (independent variable block). Instantaneous derivatives of the signaling dynamics, corresponding to the slopes of the phosphorylation dynamics signals, were also added to the X block. (3) where RSS represents the residual sum of squares of predicted Y, SSY tot.corr the total variation in the Y matrix after mean centering and scaling and PRESS, the predictive residual sum of squares defined as (4) calculated by cross validation. The appropriate number of components is defined as the optimum trade-off between goodness of fit and predictive ability. For more detailed description on PLSR modeling we refer to a previous explanation (28).  (Fig. 1b). Rapid mixing is achieved by chaotic mixing using the staggered herringbone mixers (SHM) in the channels (Fig. 1d) (20,29). We have previously showed that eight cycles of herringbones are sufficient to mix the solutes completely (20), but we observe that this configuration leads to a predictable cell focusing pattern in specific outlets, insensitive by guest on May 8, 2020 to the additional number of cycles of SHM used or the Reynolds number. To redistribute cells more evenly to facilitate the downstream cell and lysate collection, grooves with randomizer geometries and a large splitting area were added after 9 cycles of herringbones. A flow rate of 44 L/min was chosen as it led to the most uniform cell distribution in the 8 outlets (Fig. 1b).

Microfluidic chip design and operation
The second module receives the cells that have been stimulated for the desired time in the tubings. The length of the tubings was determined using the following formula: where represents the length of a tubing, t represents time spent in the tubing, Q is the volumetric flow rate and corresponds to the radius of a tubing. Time spent in the tubing is defined as (6) where (7) represents the time spent in the first module and in the pressure drop channels with rectangular cross sections (L, H and W are respectively the length, height and width of the channels) and is the total stimulation time desired.
The pressure drop channels ensure all cells are maintained at the same flow rate (Fig.   1c). The circular geometry of those channels minimizes shear forces that cells are subjected to.
The pressure drop channels did not cause measurable inertial focusing or separation of cells.
After stimulation, cells are split into two equal populations for lysis or fixation to quench the by guest on May 8, 2020 reaction. Thus, one single experiment yields 8 lysates and 8 fixed cell populations for 8 different time points, from 0.5 minutes to 7 minutes. The modularity of the chip gives the possibility to sample at various time points but since we were interested in sampling early dynamics after TCR engagement, we chose to only measure protein activation within 7 minutes of stimulation.
Each of the eight lysates was analyzed for six proteins; therefore each experiment resulted in 48 dynamic measurements. Our data suggests that flow cytometry has a larger dynamic range than our Luminex bead based assay.

Global characteristics of age-associated protein expression and activation changes in human CD8+ T cells
Although aging of immune cells has been the focus of intense research, it is difficult to reconcile reports characterizing a variety of biomarkers over a range of culture conditions. By maintaining uniform expansion and sampling conditions, we obtained the data in figure 3 presenting trends about 1,000-fold expansion, which is consistent with previous observations using this culture method (31,32). After 12 population doublings, an increase in IL-2 stimulation or in the number of beads did not allow for further growth (Fig. S1). The cell population could be maintained in culture for several weeks with appropriate culture conditions (IL-2 stimulation and fresh media), suggesting the state of replicative senescence had been reached (33). At different stages in culture, cells were sampled in order to observe changes in cell morphology, cycle, phenotypical markers and signaling as they "aged" in culture (Fig. S1).

Changes in cell morphology and cell growth arrest:
We observed a progressive decrease in cells in S/G2 phases as cells expand in culture, as others have reported (34). At the end of the proliferation phase, cells were found primarily in the G0/G1 phase of cell cycle (Fig. 4a). We also observed changes in cell size and shape as cells age through forward and side scatter detection by flow cytometry. Cell size progressively decreased with age and the variance in cellular shape increased (Fig. S2). Prior proteomic analysis of in vitro cultured CD4+ T cells identified profilin-1 as a potential biomarker for senescence (19); therefore this protein was included in our panel.
However, no statistically significant changes in profilin-1 expression were observed in our CD8+ by guest on May 8, 2020 cell population as a function of days in culture (Fig. S2). Accumulation of the cell cycle checkpoint p16 ink4 has been observed in senescent T cells (34,35). However, in our study, the increase in p16ink4 expression during in vitro aging was found to be not significant (Fig. S2).
Overall decrease in protein activation following T cell receptor stimulation: Along with changes in surface marker expression, it has been proposed that T cell function decline with age could be due to the development of defects in the transduction of mitogenic signals following T cell receptor stimulation (40). From the lysates yielded by the microfluidic chip, we quantified the levels of phosphorylated CD3, Lck, Zap70, LAT, ERK and CREB within 7 min after T cell receptor ligation using a high-throughput, multiplex bead-based assay (20). T cell receptor stimulation led to an increase in the levels of phosphorylated proteins, with the magnitude and kinetics of activation dependent on the protein, the donor and the age of the cell population ( ERK activation has been shown to be essential in mediating proliferation and telomerase activation (41,42) and displays digital or analog activation patterns depending on the strength and the nature of the stimulus (43-45); thus we investigated this protein by single cell analysis.
Analog activation appeared with anti CD3-CD28 stimulation on our CD8+ T cell population throughout the duration of culture (Fig. 2d). This analog activation is consistent with observations by Singh et al. who also relied on anti-CD3/CD28 as the means of T cell activation (44). As with the lysate dynamics, an overall decrease in signaling was also observed with time in culture. We also measured the heterogeneity in cellular ERK activation within the population, as determined by the coefficient of variation of the histograms for each time point. Higher coefficients of variations are observed very early and late with respect to culture time (Fig. S5).
Heterogeneity in the composition of the initial cell population isolated (naïve vs. effector vs. memory) might result in this high CV observed early in culture.
Extraction of biomarker combinations of cellular age in culture and prediction of cellular age and quality: Several biomarkers, such as CD28 and CD27, present modest correlation with time spent in culture; however their rate of change is very donor-specific and therefore individual markers are not sufficient to distinguish non-senescent populations with early-senescent by guest on May 8, 2020 populations. We therefore sought a combination of markers suitable for quantifying age of a cell population. Hierarchical clustering is a common method used in microarray analysis to extract clusters of genes or proteins having similar response to the same environmental factors (46). Applied to our signaling and flow cytometry data, hierarchical clustering of protein phosphorylation time courses suggested higher phosphorylation values in "young" cells (Fig. 3).
However, the clusters are very donor-specific. Combinations of biomarkers consistent for all donors do not emerge. In addition, this cluster analysis technique, very useful to interpret large dataset and group "like" variables, cannot convey quantitative contributions of markers with predictive power of T cell "age" or quality. Therefore, we applied further analytical methods for extracting markers of aging and their relative contribution to aging that would be universally predictive across donors.
A combination of dynamic signaling and expression metrics (Lck, ERK, CD28 and CD27) are the most informative markers in predicting cellular age To relate cell age to network activation and phenotypic markers, we constructed a PLSR model.
This data-driven modeling approach allows for complexity reduction in multivariate protein expression and signaling data to identify the most informative variables, and has been previously used to predict cellular fate or cytokine production (47,48 Clusters of predictor variables that are highly correlated with outcome variables can be visualized by their proximities to one another on the weight plots in principal component space (Fig. 5b). Signaling information is mostly anti-correlated to age, consistent with the general understanding of altered signaling dynamics. An analysis of the weight of each variable highlights the importance of the proteins ERK, Lck and LAT for predicting cellular age (Supplemental Table 2 phase appear to be heavily negatively weighted in the first component. Somewhat surprisingly however, the heterogeneity of the cell population in cellular shape, the heterogeneity in CD57 expression, and the basal level in phosphorylated ERK emerged as being highly positively correlated with cellular age (Fig. 5c). The heterogeneity in ERK phosphorylation in stimulated cells emerges as important in the second component (Fig. 5d).  Table 4). An analogous modeling exercise was performed with only surface marker data (Fig S6). CD27 and CD28 expression, as well as CD57 coefficient of variation emerged as the most informative metrics.
The PLSR prediction of days in culture was improved in the signaling model (R 2 Y = 0.84 vs. 0.78) but slightly less accurate in predicting population doubling (R 2 Y = 0.95 vs. 0.98).
Signaling information can predict co-stimulatory molecule expression: CD27 co-stimulation with TCR enhances cell expansion and promotes cell survival (49). The co-stimulatory molecule CD28 is involved in T cell receptor signaling amplification by PI3K activation and downstream calcium mobilization (50) as well as indirect ERK activation via Lck (51,52). Therefore, the loss of CD28 CD27 surface expression with a correlation coefficient ranging from 0.75 to 0.91 (Fig. 7c-e).
CD27 mean fluorescence intensity is highly correlated to the instant derivative of ERK phosphorylation after one minute of stimulation and anti-correlated to the slope of deactivation of CREB after 2.5 min of stimulation. Lck phosphorylation is strongly correlated to CD28 expression and the rate of decay in LAT phosphorylation after 3.5 min of stimulation is strongly anti-correlated to CD28 expression. Globally, CD28 expression clusters with early dynamics (0.5-3.5 minutes) of ERK, Lck and LAT (Fig. 7b).

Discussion
The outcome of immune based therapies, such as adoptive T cell transfer therapy or engineered vaccines, is dependent on the quality of T cell clones used (11). Progress has been achieved towards improved T cell expansion methods in the past few years (53): new culture media as well as improved stimulation techniques have been developed (54)(55)(56), and cord blood as been harnessed as a source of non-senescent lymphocytes for tumor immunotherapy (57). To assess the quality of expanded cells, researchers/clinicians generally perform functional assays (e.g. cytokine production or cytotoxic T-lymphocyte (CTL) assays (58)) or examine surface marker expression, such as the extent of CD28 loss or the appearance of late differentiation markers such as CD57 (37). Univariate, static assays measuring expression of surface markers alone are not always accurate for prediction of cell functionality. While CD28 is considered as a biomarker of immunosenescence, its expression can be downregulated by tumor-necrosis factor (TNF) (59,60) or upregulated by IL-12 in CD4+ cells (61). It was also shown that CD8+CD28-cells are able to proliferate (62) and therefore the loss of CD28 is not necessarily used statistical discriminant analysis to identify potential biomarkers (19). Each of these prior studies has relied on static information from the cells rather than functional dynamics of TCR activation. As there is no consensus for the best biomarkers of in vitro aging for CD8+ T cells, we chose to utilize a multivariate approach, combining surface phenotype of CD8+ lymphocytes and intracellular signaling.
In this work, we developed a microfluidic tool and a statistical model to evaluate cell responsiveness and accurately predict cell "age" and quality respectively. It has been observed that T cells are able to respond to stimuli within seconds, such as TCR ligation initiating a burst of calcium (69). The design of the microfluidic device enabled sampling of the rapid protein phosphorylation dynamics in the first few minutes following stimulation with minimal standard error, and allowed us to obtain accurate measurements in a high-throughput manner. As shown previously, chaotic mixing and flow in narrow channels at low Reynolds number does not elicit adverse stress response (20). The ability to reproducibly sample stimulated cells with 30-90 second intervals enabled the use of instant derivatives as additional variables in the regression model; thus, the microfluidic device provided an additional benefit in enhancing the by guest on May 8, 2020 information content from the signaling data. Although others have reported derivatives of time courses to be less informative in PLSR analysis (27,47,48,70,71), these studies relied on signaling dynamics that extended for hours rather than minutes. Because the device can stimulate a small number of cells with high temporal resolution and subsequently lyse and fix them in parallel, we could seamlessly "stitch" together complementary measurements for our statistical analysis from populations of cells treated identically, not only from day to day, but also within a particular stimulation experiment. Because donor-to-donor variability is a confounding factor in deriving robust biomarkers of senescence, the technological platform minimizes the experimental data variance so that meaningful dimension reduction could be performed, as shown by the conserved trends of signaling of phosphorylated ERK (Fig. 2). signaling in different subcellular compartments, e.g. nuclear transcription factors, not directly accessible with detergent-based lysing.
We developed multiple partial least square regression models to explore properties of the data and be able to predict physiological age of ex vivo expanded cells. PLSR models have been applied previously to understand complex signaling networks involving multiple inputs and multiple outputs, without prior knowledge of the network structure (27,47,48,70,71).
Prior application of this type of modeling approach has demonstrated that antigenic information content is encoded within downstream phosphorylation events in T cells (48). This data-driven modeling technique is particularly well suited to carry multivariate analysis and predictions in extensive datasets, and therefore could be applied to extract correlations between age, signaling, and surface phenotype expression. From the analysis, CD28 and CD27 expression emerged as key markers in determining cellular age and also correlated strongly to intracellular phosphoprotein dynamics. Interestingly, the model did not find any positive correlation between the effector memory surface markers CD57 and CD45RO and age in culture, contrary to a previous observation (37). This might be due to the culture conditions specific to each study. It has been reported that the anti-CD3 and anti-CD28 presentation as used here preferentially expand CD45RO-cells (55). CD57 is a marker associated with the endstage T cell differentiation. Elevated expression of CD57 has been observed on tumor-specific T cells. It may be the consequence of persistent chronic antigen stimulation, resulting in the accumulation of cells capable of rapidly secreting cytokines but that lack proliferative capacity (37). IL-2 supplementation may have induced an increased loss of CD28 expression (72), blocking cell proliferation before CD57 upregulation. Thus, this model is limited in scope to the by guest on May 8, 2020 culture conditions assayed yet provides a generalizable approach that could be used for other expansion methods.
Our model emphasizes the importance of early signaling dynamics (1-3.5 min after TCR engagement) to explain and predict cellular age. Others have reported alterations in early signaling events in old mice (73) and in human T lymphocytes from elderly patients (74), furthering the common characteristics between ex vivo culture expansion and immuosenescence. Three major proteins, Lck, LAT and ERK, were extracted from the model as having primary roles in accurately predicting cellular age. Lck is the first protein to be phosphorylated after TCR engagement, leading to phosphorylation of the adaptor protein LAT via Zap70, and downstream phosphorylation of ERK. Impaired redox regulation (75), reduced calcium release (73,76) and altered membrane rafts composition (77) are possible explanations for the impaired activation of those proteins.
Previous studies have reported higher cell-to-cell variation with increasing age (78,79); however the modeling results suggest that this trend is more nuanced and depends on the protein in question. We find that the heterogeneity in CD57 protein expression and cellular shape is correlated with cellular age. In contrast, no direct correlation between age and cell-tocell variation of ERK phosphorylation was discovered (Fig. S5), and yet as part of a multivariate analysis this CV variable still possesses a high predictive power. This suggests robustness in ERK response to intracellular noise created with aging (80). The multivariate PLSR model is able to assess the age and therefore the quality of a cell population for cells expanded using CD3-CD28 bead stimulation with IL-2 supplementation. Past by guest on May 8, 2020 multivariate models associated signaling information with cellular fate (70), cytokine production (48), or drug response (81). These models enabled the creation of new hypotheses to test the validity of the model. In this application, aging is not a cell fate that we can modify and the biomarkers extracted are not correlated in a causal fashion to the age of the cells. This lack of a testable mechanism poses a limitation to the study; however, the model yields novel insight on the most informative markers of the array selected for sampling, and predicts cellular age from those specific markers. We also envision this model to be a possible diagnostic tool to quantify immune age of elderly individuals or individuals presenting accelerated immunosenescence, such as HIV patients.
In summary, the design of a novel microfluidic chip enabled the statistical analysis of T cell aging by minimizing error in the sample handling between days and across multiple donors.
We took advantage of chaotic mixing geometries and a modular chip design to capture the early phosphorylation events associated with TCR ligation, which in turn proved to be highly informative in the partial least square regression prediction of population doubling and days spent in culture. Our findings point to a cell signaling-based assessment method that could quickly evaluate patient-derived cells for degree of population doubling. The general approach described here, combining fixation and lysing on-chip, can facilitate the integration of single-cell information with population-averaged techniques such as multiplexed immunoassays or mass spectrometry.
Acknowledgements by guest on May 8, 2020   Predi g Cytotoxic T cell Age Figure 3: Complete dataset. For each donor, dendrograms were generated with signaling measurements from lysates (le and middle), protein phosphoryla on qua ca and instant deriva ves of those measurements. Flow cytometry dendrograms (right) contain ERK phosphoryla on, cell cycle, cell morphology, surface marker, profilin-1 and p16 ink4 expression. Each row of 3 dendograms corresponds to a different donor. In each dendogram, a row corresponds to a sampled me. Each column corresponds to a par cular variable measured.