Mathematical model of a personalized neoantigen cancer vaccine and the human immune system

Cancer vaccines are an important component of the cancer immunotherapy toolkit enhancing immune response to malignant cells by activating CD4+ and CD8+ T cells. Multiple successful clinical applications of cancer vaccines have shown good safety and efficacy. Despite the notable progress, significant challenges remain in obtaining consistent immune responses across heterogeneous patient populations, as well as various cancers. We present a mechanistic mathematical model describing key interactions of a personalized neoantigen cancer vaccine with an individual patient’s immune system. Specifically, the model considers the vaccine concentration of tumor-specific antigen peptides and adjuvant, the patient’s major histocompatibility complexes I and II copy numbers, tumor size, T cells, and antigen presenting cells. We parametrized the model using patient-specific data from a clinical study in which individualized cancer vaccines were used to treat six melanoma patients. Model simulations predicted both immune responses, represented by T cell counts, to the vaccine as well as clinical outcome (determined as change of tumor size). This model, although complex, can be used to describe, simulate, and predict the behavior of the human immune system to a personalized cancer vaccine.


B.1. Estimating peptide concentration in a vaccine from amino acid sequences
We extracted the amino acid sequence of the peptides used for each patient's personalized vaccine used in the melanoma study [1] from their supplementary material. First, using all amino acid sequences of peptides in each vaccine, we calculated the protein molecular weight (KDa) by using the online tool The Sequence Manipulation Suite: Protein Molecular Weight [2]. Then, in [1], it was stated that each peptide had a weight of 0.3 mg. We used the online tool [3] to covert the weight concentration to pmol. In Table A1, we lay out each of the components necessary to calculate the peptide concentration, Dose p in pmol for input in our model. Table A1. Peptide dose conversion from mg to pmol by patient with melanoma in [1].

B.2. Dendritic cells
Carrying capacity, K DC . To estimate the carrying capacity of immature dendritic cells in the total volume of the subcutaneous tissue injection site, we used the total number of dendritic cells in Figure 4A of [4]. Since this article shows the total number of DCs by gender, we use the total number of DCs in females (approx. 27 number of cells/µL) as the lower bound of our range and we use the total number of DCs in males (approx. 35  (1) below) of the second term in the equation of our model, D I (t), representing the maturation rate of immature DCs to data from Figure 6A in [5] using Mathematica 12.0 [6] by looking for the best fit. The data from Figure 6A in [5] was extracted using WebPlotDigitizer [7]. Figure 6A in [5], shows how DCs mature as a function of poly(I:C) adjuvant. In this figure, authors use CD83 since it is known that these are activation markers for antigen presenting cells [8]. We re-scaled %CD83 + cells to proportion of matured DCs in the y-axis and microgram per milliliter to milligram per liter in the x-axis of Figure 6A resulting in Fig 1 (see below).
Using separation of variables method, we find the analytical solution required to estimate r D and K a , and integrate such that for all t ∈ [0, T ), According to the experiments done in [5], the cells were cultured for 24 hours before taking measurements of maturation. Hence, we set T = 1 (1 day) and the function to be fitted is: 3/12 Using the built-in function 'NonlinearModelFit' in Mathematica [6], we obtain the following parameter estimates from equation 1:  Figure 6A in [5]. 90% confidence bands for mean predictions are in blue. [90% CI]* denotes modified CI by replacing negative lower endpoint with zero as described in [9].

B.3. Estimation of parameters corresponding to the model equations that represent naïve T-cell populations
Initial counts of naïve T-cells. In Ott et al. [1], eligible patients needed to have a lymphocyte count of at least 800 cells per microliter. We use these number as a baseline to estimate naïve T-cell counts at t = 0 in our simulations. The average adult weighting 58 to 80 kilograms has about 4.5-5.7 L of blood. If we assume that patients have an average of 5 L of blood, then each patient has at least 4 × 10 9 lymphocytes. On average, using standard methods such as Ficoll-Paque density gradient centrifugation for PBMC yields 0.5-2×10 9 cells/L of blood [10]. This will give us a range of (2.5, 10) × 10 9 cells, with a median of 6.25 × 10 9 cells.
In Bittersohl et al. [10], authors provide percentages of different cells types which compose the PBMCs in a healthy adult human as shown in Table 9.1 of [10] (see Table A2 below). From this table, we obtain that 70% of PBMCs are lymphocytes for which 60 out of 70 percent (about 86% of total lymphocytes) are T-cells and 10 out of the 70 percent (about 14% of total lymphocytes) are B-cells. With these percentages, we obtain a T-cell range of (2.15, 8.6) × 10 9 with a median of 5.38 × 10 9 T-cells. Since this median is larger than the baseline calculated above, we assign 5.38 × 10 9 cells as the initial naïve T-cell count in our simulations.
T-cell Carrying capacity. In order to have a one-to-one comparison of T-cell count from extracted data and simulation, we first estimate the average number of PBMCs in a human. Following [11], we assumed that a human has 1.22 × 10 12 lymphocytes. Using Table A2, we estimate there are 1.43 × 10 12 PBMCs. Thus, we conclude that T-cells = 8.57 × 10 11 (60% of PBMCs).

B.4. IFN-γ ELISPOT of PBMCs data extraction
In Ott el al. Figure 2b [1], the authors provided longitudinal plots for ex-vivo assay responses to peptide pools for each patient in the trial in units of per million PBMC. These data points were obtained by using an IFN-γ ELISPOT assay. This method allows the quantification of the number of CD4 + or CD8 + T-cells which secrete IFN-γ in response to a stimulation with a specific antigen or peptide [12,13]. Moreover, each spot within an ELISPOT well identifies a single cell that has released a measurable amount of cytokine, i.e., the number of spots is a direct measure of the frequency of cytokine-producing T-cells in the PBMC population [12]. Thus, it is a reasonable assumption that the number of spots to number of cells is a one-to-one conversion.
To extract data from plots and images, we used WebPlotDigitizer [7] and plotted the extracted data (see Fig 2A) by patient. Each patient's frame contains the data of the four pools (Pool A -Pool D) which make up one vaccine dose and a baseline data (mock) . Fig 2A is plotted in the form of cell counts by multiplying each extracted data point by 1.43 × 10 6 (see section: T-cell Carrying capacity), since each data point is in spot-forming cells per 10 6 PBMCs. Then, to plot the total cell count of ex-vivo assay responses to peptides vaccine for each patient, we summed the data points of Pool A through D each subtracted by the Mock data at each time point and arrive at the numbers shown in Fig 2B. Due to the stochastic nature of ex-vivo assays, one caveat to subtracting the Mock data is that it may lead to negative cell counts on the treatment initiation day for which we assigned a small (10 −6 ) number of activated cells. Additionally, we assumed there is a +/-15% variation to the ELISPOT counts presented by Ott et al., as it has been shown in [14] that the relative experimental error for 0-200 number of spots counted per well is between 0.1-0.2. CD4 + and CD8 + T-cell initial counts. The first data point at day 0 for each patient was used to assign to each patient's initial CD4 + and CD8 + counts. Given that these data only tell us the T-cell count without making distinction between cell subtypes, based on Table A2, we assume that T-cell counts on day 0 are distributed as 70% CD4 + T-cells and 30% CD8 + T-cells.

B.5. Parameters corresponding to the equation describing tumor cell population
Estimating maximum cancer growth rate, r. In Liu et al. [15], authors studied the rate of growth in melanomas, finding that one third of the melanomas grew at least 0.5 mm per month. The median monthly growth rate for nodular melanomas was 0.49 mm in diameter, superficial spreading melanomas at 0.12 mm in diameter, and lentigo melanomas at 0.13 mm/month [15]. The patients sample for our simulations in [1] were diagnosed at different stages of melanoma cancer (see Table A4). Four of the patients who did not have recurrence after vaccine initiation had nodular melanomas (i.e., the tumor has spread to nearby lymph nodes [16]). However, all patients before initiating vaccine treatment, underwent surgery, thus, we assumed the growth rate of any residual after surgery is 0.15 mm/month in diameter. That is, we are assuming the best case scenario for these patients, regardless their initial tumor stage diagnosis.
To determine the number of cancer cells of a tumor in a spherical shape, we considered that each cell is about 20 µm in diameter and a 1-mm cancer has about 100 thousand cells [17]. Therefore, to directly calculate the cell count in a tumor of diameter d, we used the formula of the volume of a sphere modified to incorporate a 25.95% of void volume [18], and multiply by where d is the diameter in millimeters (mm).
The tumor growth rate is then calculated as follows: Estimating initial tumor cell count, T (0), by patient. To estimate the initial tumor cell count on day 0 (t = 0), we used the information provided for each of the patients in the 'Extended Data Table 1: Characteristics of patients' along with the clinical event timeline from surgery until time of data cut off in ' Figure 1', both in [1]. Then, following the guidelines from the American Cancer Society [16] and Melanoma Research Alliance [19] for melanoma skin cancer staging and using nodal staging and metastasis information on [20] and [21], respectively, we made the following assumptions which are then summarized in Table A3: (a) Given the tumor diameter ranges in [19] for each of the tumor thickness (T) in column 2 of Table A3, we make a definitive assumption for the tumor diameter.
(b) According to [20], lymph nodes are considered malignant if they measure more than 1 cm in the short axis diameter. Though, the size threshold varies with anatomic site and underlying tumour type. For example, in rectal cancer lymph nodes are considered pathological when they measure more than 5 mm. Hence we assigned (column 5 of Table  A3) the number of tumor-involved nodes according to the code (N).
(c) If a patient has metastasis (M), then in [21] it was stated that lung metastasis ranges between 3-7 cm. Hence, we assumed that patients with metastasis will have three 5-cm in diameter lung metastasis (column 7).
As an example, we will walk you through the calculation to obtain the tumor cell count of Patient 2 at the beginning of cancer treatment (column 6 in Table A4):  Table A3. Definition of tumor diameter according to the American Joint Committee on Cancer (AJCC) TNM system [20].
• Patient 2 in [1] was diagnosed with T4N0M1b. Using Table A3, we conclude that T4 corresponds to a tumor diameter of 5 mm; N0 means there are no tumor-involved nodes; M1b corresponds to one 50-mm lung metastasis. Then using equation (2), we compute the tumor cell count for a tumor diameter of 5 mm, adding three times the tumor cell count for a tumor of 50 mm in diameter, Then, the tumor cell count after surgery is considered to be 5% of the tumor size before surgery. Thus, we multiply column 3 by 0.05. Lastly, we use the equation for tumor cells (Eq. (15) in the main text) without the second term which represents the intervention of the immune system. This means, we assume the tumor cells continue to grow without any mediation from the immune system.
In the equation above, we use 0.004 as the maximum growth rate estimated previously and 1.45 × 10 10 as the tumor carrying capacity since this is the largest tumor cell count before surgery among patients (see column 2 in Table A4 Table A4. Estimation of initial tumor cell count by patient at the beginning of vaccine treatment. Third column is estimated using Table A. Fourth column is estimated using equation (2). Fifth column is the time lag between surgery and initiation of cancer vaccine treatment. Sixth column is estimated using equation (3).

B.7. Detailed methodology for global sensitivity analysis
We chose N = 1000 as the number of LHS/PRC simulations. Following the approach and algorithm described in [23], we constructed an LHS matrix of size N by K (N : no. of rows; K: no. of columns). Hence, each row corresponds to one simulation and contains a randomly selected and equally distributed value for each of the uncertain parameters of interest. Each of these rows in the LHS matrix is then used to calculate each of the outcome variables of interest producing N observations that are then used to assess the sensitivity of the outcome variables to the uncertainty in the input parameters. Next, we calculate the PRCC value for each outcome variable at specific time-points of interest (t = 24, 111, 146) following the approach of [23]. Sensitivity analysis simulations were performed in Mathematica, Version 12.0 [6].

B.8. Additional uncertainty and sensitivity analysis for model outputs
We briefly present the results of the global sensitivity analysis performed for two model variables of interest, activated T-cells and tumor cells, using parameters found in the literature which may have associated uncertainty that could affect model outcome, such as error in measurements or natural variations. As observed in Fig 3, it was found that both model variables of interest, activated T-cells and tumor cells, are highly sensitive to parameters α p , Λ and δ M (internalization rate of peptides by DCs, maximum growth and death rate of mature DCs, respectively). The sensitivity produced by these parameters are shown to decrease throughout time. Individually, in addition to α p , Λ and δ M , model variable for activated T-cell was found to be highly sensitive to b 4 , σ 4 , µ 4 and µ 8 (maximum growth of naïve CD4 + T-cells, maximum activation rate of CD4 + T-cells, and death rate of activated CD4 + and CD8 + T-cells, respectively). The model variable for tumor cells, was found to be highly sensitive b 8 , σ 8 and ρ 8 (maximum growth of naïve CD8 + T-cells, maximum activation rate CD8 + T-cells, and proliferation rate for activated CD8 + T-cells).

B.9. Methods for the 'Application' section
Using all patient-specific parameters ( Table 2 and parameters obtained from [1]), we created six different patient profiles, and varied their initial tumor cell counts, which are assumed to be contained in a spherical shaped tumor with diameter d (see Table A5). To create a wide range of initial tumor cell counts corresponding to different tumor sizes (sample size: 2000), we used the LHS technique to get a uniform and equal distribution of tumor cell counts ranging from 30,000 to 2 × 10 10 cells. Notice that this set of different initial conditions implicitly correspond to a total number of cells x number of days (or months) after having resection surgery. Keep in mind that there may be at least 5% tumor residue and may continue to grow post-surgery. Therefore, the estimated tumor cell count (or diameter) in Table A5 is the tumor cell count at day 0 of vaccine treatment.