Model fitting for pERK data analysis

For the single-cell analysis, flow cytometry measurements of human PBMC were first gated on CD8⁺ PD-1⁺ cells using the R package flowCore (Hahne et al., 2009). Quantitative single-cell fluorescence levels of pERK, CD28, and PD-1 were then log-transformed and standardized. A linear Bayesian statistical model was fit to determine how quantitative CD28 (x₁) and PD-1 (x₂) expression impacts phosphorylation of ERK (y) 5 min after stimulation with APC that either do or do not express PD-L1 (x₃). A quadratic term for PD-1 expression was also included to account for mild nonlinearities in its effect. This model decisively outcompeted an analogous model without the quadratic term according to the widely applicable information criterion . Thus, the final model was as follows:

This model was fit using the R package brms (Bürkner, 2017) using 5,000 Markov chain Monte Carlo iterations, including 1,000 warm-up iterations. Weakly regularizing Gaussian priors were centered at 1 for β₀ and 0 for all other parameters. After model fitting, 10,000 samples were drawn from the model’s posterior distribution and used to predict the expectation and 95% credible interval for pERK expression as it depends on CD28, PD-1, and PD-L1 expressions.