A Tumor-Immune Interaction Model for Synergistic Combinations of Anti PD-L1 and Ionizing Irradiation Treatment

Combination therapy with immune checkpoint blockade and ionizing irradiation therapy (IR) generates a synergistic effect to inhibit tumor growth better than either therapy does alone. We modeled the tumor-immune interactions occurring during combined IT and IR based on the published data from Deng et al. The mathematical model considered programmed cell death protein 1 and programmed death ligand 1, to quantify data fitting and global sensitivity of critical parameters. Fitting of data from control, IR and IT samples was conducted to verify the synergistic effect of a combination therapy consisting of IR and IT. Our approach using the model showed that an increase in the expression level of PD-1 and PD-L1 was proportional to tumor growth before therapy, but not after initiating therapy. The high expression level of PD-L1 in T cells may inhibit IT efficacy. After combination therapy begins, the tumor size was also influenced by the ratio of PD-1 to PD-L1. These results highlight that the ratio of PD-1 to PD-L1 in T cells could be considered in combination therapy.


Introduction
Immunotherapy using immune checkpoint blockade (IT) is an anti-cancer therapy that recovers immunity by suppressing various tumor mechanisms that evade the immune response [1][2][3][4]. In particular, new ITs primarily aim to inhibit cytotoxic T-lymphocyte-associated protein 4 (CTLA-4), programmed cell death protein 1 (PD-1), and programmed death ligand 1 (PD-L1) to enable T lymphocytes to attack cancer cells [5][6][7]. Despite the impressive potential of these checkpoint blockade in the treatment of various cancers, this therapy strategy remains challenging because the response to a certain subgroup of IT is varied among the cancer patients [8,9]. In recent years, several studies have aimed to achieve synergistic effects through combination therapy to compensate for the shortcomings of monotherapy [10][11][12][13].
There are mathematical models of immune-tumor interactions to demonstrate the data of Deng et al. [14]. These studies [15,16] revealed the synergistic effect of anti-PD-L1 and IR combination therapy in mice. Chappell et al. [15] discussed that IT is directly involved in increasing T cell levels and indirectly stimulates tumor death via T cells. The underlying assumption considered that IR affects the Furthermore, 200 μg IT is administered every three days for a total of four doses in a 50 μL volume after beginning IT on day 14.
The second term of the right-hand side (RHS) indicates tumor elimination by IT and IR therapy.
The first term of the RHS indicates that the PD-1-PD-L1 complex suppresses T cell growth. The second term of the RHS describes natural cell death and additional elimination by IR.
The first term of the RHS represents the increase in PD-L1 expression in tumors and T cells. The second term indicates inhibition of PD-L1 by IT. The third and fourth terms indicate the binding process of PD-1 and PD-L1. We assume that this reaction occurs according to the mass action law, indicating a one to one correspondence.
In the first term of the RHS, , 0 ≤ ≤ 1 determines the ratio of PD-1 to PD-L1 induced by T

Data Derivation from the Study
Published data was obtained from the experiments of Deng et al. Their study was to verify the synergistic antitumor effects of IR and IT. A combination of IT and IR significantly enhanced the inhibition of TUBO growth in their study. The outline of the experiment process was as follows. BALB/c mice were inoculated by subcutaneous injection with 1 × 10 6 TUBO cells (the size 5 mm 3 ). After injected s.c. and waiting for 14 days (growing), tumors were locally treated with single 12 Gy of IR and 200 µg IT (antibody, clone 10F.9G2) by intraperitoneal injection every three days for a total of four doses. Tumor size was measured in mice with no therapy, IT only, IR only, and combination therapy with IT and IR. For reproducing the experimental data, the data values were extracted from the study. Adobe Illustrator is used to find the coordinates of published data. Data values, shown in Figure 2 (left top panel), are presented in Table 1. Deng et al.'s study was to verify the synergistic antitumor effects of IR and IT. A combination of IT and IR significantly enhanced the inhibition of TUBO growth in their study. conducted with the nonlinear least square method using MATLAB and Berkeley Madonna. Using the estimated parameters, the model captured the synergistic effect of the combination therapy consisting of IR and IT, as shown in Figure 2 (IR + IT). Additionally, the model failed to capture data around day 22. This is because repeated combination therapies cause a slightly overstated discrepancy in the model. However, the overall data fitting is reliable, and the model captures the efficacy of the combination therapy in mice. In the control case, tumor size, PD-1 and PD-L1 expression were positively exponentially proportional. In this case, the regression (function) over time from day 0 to 35 was ( ) = 11.45 . ⋅ with = 0.9975. Cases were compared before and after therapy. Before therapy, the

Mathematical Modeling
Let C, T, PL, PD, S, R, and A denote tumor (volume), the density of CD8+ T cells (T cells), PD-L1, PD-1, PD-1-PD-L1 complex (µg/µL), IR (Gy), and IT (µg/µL), respectively. The tumor grows logistically, and T cells are activated by the cytokines interleukin 12 (IL12) and are proliferated by interleukin 2 (IL2). Both IL12 and IL2 are activated by stimulation of the tumor. Among the several types of immune cells playing a role against cancer, we have focused only on T cells. IL-12 and IL-2 play a significant role in the activation and differentiation of the population of T cells [17]. The tumor expresses PD-L1 on the cell surface, and T cells express both PD-L1 and PD-1 on their surface. The binding process of PD-L1 and PD-1 induces negative feedback in T cells. IT inhibits the growth of both tumor and T cells and inhibits PD-L1 by preventing the formation of the PD-L1-PD-1 complex, which indirectly prevents the inhibition of T cell growth.

Modeling Therapy with IR (R) and IT (A)
An exponential decay model was used to model for IR as follows: In the case of IT, a pharmacokinetic model was used with a two-compartment model as follows: Pharmaceutics 2020, 12, 830

of 17
Furthermore, 200 µg IT is administered every three days for a total of four doses in a 50 µL volume after beginning IT on day 14.

Compartmental
The second term of the right-hand side (RHS) indicates tumor elimination by IT and IR therapy.
The first term of the RHS indicates that the PD-1-PD-L1 complex suppresses T cell growth. The second term of the RHS describes natural cell death and additional elimination by IR.
The first term of the RHS represents the increase in PD-L1 expression in tumors and T cells. The second term indicates inhibition of PD-L1 by IT. The third and fourth terms indicate the binding process of PD-1 and PD-L1. We assume that this reaction occurs according to the mass action law, indicating a one to one correspondence.
In the first term of the RHS, η, 0 ≤ η ≤ 1 determines the ratio of PD-1 to PD-L1 induced by T cells. That is, PD-L1 is more expressed as η approaches one.

Model Reduction: Quasi Steady-State Approximation (QSSA)
We assume that PD-1-PD-L1 complex S is quickly saturated. Then equilibrium is reached within a short period of time, and so we may apply quasi-steady-state assumption (QSSA) and obtain k on PL · P = k o f f S. Assume that k D = k on k o f f . Then the system is simplified as follows.
Activation by IL12 Growth by tumor and T cells Pharmaceutics 2020, 12, 830 Growth and the ratio expression level of the T cells Descriptions or units of the estimated parameters and compartments are presented in Table 2. Table 2. Estimated parameters and initial values used in the model are presented together with their descriptions and units.

Initial/Parameters Description and Units Estimated Values
Initial lymphocytic density of CD8+ T cells (10 9 cells/L) Irradiation ( Delay rate with the mean duration 1/k 1 (1/day) 0.15 k R Elimination rate of ionizing irradiation (1/day) 0.09

Simulation using the Cancer-Immune Model with IR and IT Therapy
Combination therapy with IT and IR significantly enhanced the inhibition of tumor growth in mice, as shown in Figure 2 (top left). The mathematical model consisting of (1)-(6) was used to fit the data, as shown in Figure 2 (top right). Data from the control (no therapy), IT only, and IR only samples were fitted. From three data fittings, we estimated parameter values, as shown in Table 2. MATLAB, MathWorks ® , was used for model implementation with ODE45, and parameter estimation was conducted with the nonlinear least square method using MATLAB and Berkeley Madonna. Using the estimated parameters, the model captured the synergistic effect of the combination therapy consisting of IR and IT, as shown in Figure 2 (IR + IT). Additionally, the model failed to capture data around day 22. This is because repeated combination therapies cause a slightly overstated discrepancy in the model. However, the overall data fitting is reliable, and the model captures the efficacy of the combination therapy in mice.
In the control case, tumor size, PD-1 and PD-L1 expression were positively exponentially proportional. In this case, the regression (function) over time from day 0 to 35 was f (t) = 11.45e 0.0052·t with R 2 = 0.9975. Cases were compared before and after therapy. Before therapy, the linear regression was f (t) = 0.267t + 22.27 with R 2 = 0.8019. Linearity was well followed after therapy begins (day 14) and the regression was f (t) = 2.461t − 859.8 with R 2 = 0.8732. In the case of IT, there was a positive correlation after initiating therapy. The linear regression was f (t) = 1.451t − 315.7 with R 2 = 0.9891. Although the growth rate was suppressed by IT, the relationship between IT and tumor growth still had a positive correlation. Likewise, in the case of IR, the growth rate was lower, and the linear regression was f (t) = 1.294t − 343.4 with R 2 = 0.6277, which shows positive correlation. However, there was a negative correlation after beginning combination therapy, in which the tumor size decreases due to the synergistic effect. The linear regression was f (t) = −0.2029t + 150.7 with R 2 = 0.9731.

The Expression Levels of PD-1 and PD-L1
We utilized the model to explore which parameters cause tumor suppression among PD-1 and PD-L1. That is, we analyzed the changes in tumor size vs. parameters related to the changes in PD-1 and PD-L1. From (5)-(6), expression levels ν and µ determine the magnitude of the expression levels of PD-1 and PD-L1 on the surface of tumors and T cells. Given that T cells can express both PD-1 and PD-L1, the ratio η, 0 ≤ η ≤ 1, determines the ratio of the expression levels of PD-1 and PD-L1 in T cells. That is, if η is increasing, then so is the expression level of PD-L1 in T cells. Analyzing these parameters enables us to evaluate whether the concentration is more influential in removing the tumor among the expression levels of PD-1 and PD-L1. Specifically, we investigated how the IT efficacy is maximized depending on the expression level ratio of PD-1 to PD-L1. For this analysis, we attempted local changes in parameters, as shown in Figure 3, and two of the three parameters were compared. When modulating µ and ν (left panel), an increase in u and ν stimulates tumor growth. This is because an increase in PD-1 and PD-L1 under a fixed η causes T cell suppression, resulting in tumor growth. When modulating ν and η (middle), tumor growth depends on an increase in ν and decrease in η. This change indicates that IT could be more efficient when PD-1 is more highly expressed than PD-L1. When η and µ are varied (right), tumor growth depends on decreasing η and increasing µ. Thus, in the local sense, the decrease in η and the increase in ν and µ result in tumor growth. This indicates that IT efficacy is more potent when PD-1 level on T cells is smaller.
In Figure 4, we investigated tumor size vs. the sum of PD-1 and PD-L1 expression. We noticed that the total expression level of PD-1 and PD-L1 does not positively correlate tumor size when the tumor is inhibited by combination therapy, unlike IT and IR only. We inferred that this was because the balance of PD-1 and PD-L1 was broken, leading to a reduced concentration of the complex formed by PD-1 binding to PD-L1 and ultimately blocking T cell inhibition. However, these local changes in parameters could have biased influences given that these processes naturally assume that the other parameters remain constant. In particular, the change in expression levels of PD-1 and PD-L1 are simultaneously affected by various parameters.   Performing a two-parameter comparison could induce prediction bias. Thus, global sensitivity analysis is required to analyze how all parameters associated with the changes in PD-1 and PD-L1 are simultaneously changed.
From Figures 3 and 4, the model verified that the magnitude of tumor reduction was changed according to the change in the parameter, and the change in tumor size due to the change in the sum of PD-1 and PD-L1 expression was compared. The model showed that the concentration of the sum of PD-1 and PD-L1 was not proportional to tumor size after therapy begins, indicating that if tumor size decreases, the correlation between them is deregulated. We believe that the synergistic effect of combination therapy accelerated the dysregulation of this correlation. However, further studies on Performing a two-parameter comparison could induce prediction bias. Thus, global sensitivity analysis is required to analyze how all parameters associated with the changes in PD-1 and PD-L1 are simultaneously changed.
From Figures 3 and 4, the model verified that the magnitude of tumor reduction was changed according to the change in the parameter, and the change in tumor size due to the change in the sum of PD-1 and PD-L1 expression was compared. The model showed that the concentration of the sum of PD-1 and PD-L1 was not proportional to tumor size after therapy begins, indicating that if tumor size decreases, the correlation between them is deregulated. We believe that the synergistic effect of combination therapy accelerated the dysregulation of this correlation. However, further studies on the correlation between tumor size and the sum of PD-L1 and PD-1 expression should be conducted with more data.
The challenge is to determine how uncertainty in the output of a model can be apportioned to different sources of uncertainty in the model input (factors) [18]. Local methods analyze sensitivity around some (often optimal) points in the factor space, whereas global methods attempt to analyze variability across the full factor space. Changing one factor at a time means that whatever effect is observed on the output is due solely to that factor. Conversely, a non-zero effect implies influence; i.e., it never detects uninfluential factors as relevant. Local sensitivity has the above disadvantage; therefore, we performed global sensitivity analysis and analyzed tumor changes under the influence of other variables. The elimination rate k el is contained to evaluate the influence. Providing k el enables comparison to the influence of a dummy-parameter not directly related to PD-1 and PD-L1.

Global Sensitivity Analysis of the Parameter Space: Latin Hypercube Sampling Partial Rank Correlation Coefficient
We evaluated the influence of the parameters on the tumor size using global sensitivity analysis. Influences of the parameters were measured at days 10, 14, 20, and 35, using the partial rank correlation coefficient (PRCC) [18,19] based on Latin hypercube sampling (LHS) [20,21]. LHS is a sampling method that provides an unbiased estimation of the output of a model while requiring fewer samples to achieve the same accuracy of simple random sampling. For parameters x j , j = 1, 2, · · · , k and tumor y, a correlation coefficient (CC) between x j and y is calculated as where Cov(x, y) is the covariance between x j and y, and var x j is the variance of x j . CC varies between −1 and 1. Pearson correlation coefficient (PCC) characterizes the linear relationship considering residuals x j −x j and y −ŷ, with x j and y being first-rank transformed parameters and tumor size is determined by LHS, respectively. PRCC is a robust sensitivity measure for monotonic relationships between input and output. PRCC performs a partial correlation analysis on rank-transformed data in two steps: (i) input and output data are first rank-transformed, and then (ii) the linear regression model is applied. LHS is used for parameter sampling, in which each parameter is assumed to be uniformly distributed. We explored four parameters related to the changes in PD-1 and PD-L1 stimulated by the tumor and T cells. Herein, tumor size over the parameters is plotted using LHS, as shown in Figure 5. By varying the parameters associated with the changes in PD-1 and PD-L1, tumor size increases or decreases. The number of sampling parameters was 1000, and the ranges of the parameters were as follows: where U is uniform distribution. The resultant PRCCs with p-values are measured in Figures 6-9 using LHS. This was done at the following timepoints: day 10 (before therapy), 14 (beginning of therapy), 20 (during therapy), and 35 (after therapy). Additionally, PCC was calculated for comparison. PCC is intuitive because there is no change in parameters and tumor size, but the linear relationship of PRCC between variations in parameters and tumor size becomes more apparent compared to PCC.
Pharmaceutics 2020, 12, 11 of 17 Figure 5. Tumor growth curve with combination therapy of IR and IT is plotted with the parameter variations using Latin Hypercube Sampling. One thousand cases with four parameters randomly extracted from (7) were conducted. By varying these parameters, some tumors are eliminated, some reach maximal tumor sizes, and the others are in between the two extremes. The amount of combination therapy and other parameter values were unchanged except for , , , and .  Tumor growth curve with combination therapy of IR and IT is plotted with the parameter variations using Latin Hypercube Sampling. One thousand cases with four parameters randomly extracted from (7) were conducted. By varying these parameters, some tumors are eliminated, some reach maximal tumor sizes, and the others are in between the two extremes. The amount of combination therapy and other parameter values were unchanged except for η, ν, µ, and k el .
Pharmaceutics 2020, 12, 11 of 17 Figure 5. Tumor growth curve with combination therapy of IR and IT is plotted with the parameter variations using Latin Hypercube Sampling. One thousand cases with four parameters randomly extracted from (7) were conducted. By varying these parameters, some tumors are eliminated, some reach maximal tumor sizes, and the others are in between the two extremes. The amount of combination therapy and other parameter values were unchanged except for , , , and .

Analysis of Changes in Parameters
We compare PCC and PRCC with p-values, as shown in Figures 6-9. PCC directly compares tumor size vs. parameters. While PCC is intuitive, it has a disadvantage in that the relationship between parameters and tumor size is not clearly identified. Therefore, PRCC is performed to examine the influence of parameters more clearly by transforming the value. Ultimately, the global sensitivity analysis investigates the changes in essential outputs caused by simultaneous changes in parameters. We explored the parameter influence at each critical timepoint. In particular, the change in parameters can be compared before and after therapy. The influence of and did not significantly affect both tumor growth and inhibition at day 10, before therapy. In the case of , negative linearity is observed in PRCC as compared to that in PCC. In the case of , the linearity is not observed in PCC, but PRCC captures positive linearity. When combination therapy begins on day 14 (Figure 7), , and are not correlated to the change in tumor size. Tumor size has a positive linear relationship to , unlike before therapy. In Figure 8, which shows day 20 of IT therapy, tumor size and have a positive relationship, which becomes more apparent compared to PCC. After therapy (day 35), and are positively correlated with tumor size, as shown in Figure 9. This correlation indicates that the elimination rate of IT is influential to the change in tumor size when the amount of IT is small. From global sensitivity analysis, we may conclude that parameter influence over tumor size is different before and after therapy. These levels accelerated tumor inhibition after therapy began. Particularly, determining the ratio of PD-1 to PD-L1 does not follow the linear relationship to tumor size.
Changes in tumor size with respect to ν before and after therapy were investigated. There was a negative correlation in the absence of therapy (before day 14), as shown in Figures 6-9. However, a positive correlation was observed after IT, IR, and combination therapy began. This seemingly contradictory result demonstrated that, prior to initiating therapy, PCC showed a small change in tumor size (less than 3E-5), so the negative correlation may not have a significant effect on the change in tumor size. PRCC showed this difference to be more exaggerated, indicating that it was essential to examine both PCC and PRCC. Changes in after 14 days of therapy were positively correlated with tumor size. This relationship was due to the negative feedback by which the increase in caused by T cells leads to an increase in PD-L1, and a decrease in T cells by increasing the amount of PD-1 and PD-L1 binding. This mechanism increased tumor size and suggested that high PD-L1 expression levels in T cells could inhibit the efficiency of IT. Likewise, an increase in expression level

Analysis of Changes in Parameters
We compare PCC and PRCC with p-values, as shown in Figures 6-9. PCC directly compares tumor size vs. parameters. While PCC is intuitive, it has a disadvantage in that the relationship between parameters and tumor size is not clearly identified. Therefore, PRCC is performed to examine the influence of parameters more clearly by transforming the value. Ultimately, the global sensitivity analysis investigates the changes in essential outputs caused by simultaneous changes in parameters. We explored the parameter influence at each critical timepoint. In particular, the change in parameters can be compared before and after therapy. The influence of k el and ν did not significantly affect both tumor growth and inhibition at day 10, before therapy. In the case of ν, negative linearity is observed in PRCC as compared to that in PCC. In the case of η, the linearity is not observed in PCC, but PRCC captures positive linearity. When combination therapy begins on day 14 (Figure 7), k el , η and µ are not correlated to the change in tumor size. Tumor size has a positive linear relationship to ν, unlike before therapy. In Figure 8, which shows day 20 of IT therapy, tumor size and ν have a positive relationship, which becomes more apparent compared to PCC. After therapy (day 35), ν and k el are positively correlated with tumor size, as shown in Figure 9. This correlation indicates that the elimination rate of IT is influential to the change in tumor size when the amount of IT is small. From global sensitivity analysis, we may conclude that parameter influence over tumor size is different before and after therapy. These levels accelerated tumor inhibition after therapy began. Particularly, η determining the ratio of PD-1 to PD-L1 does not follow the linear relationship to tumor size.
Changes in tumor size with respect to ν before and after therapy were investigated. There was a negative correlation in the absence of therapy (before day 14), as shown in Figures 6-9. However, a positive correlation was observed after IT, IR, and combination therapy began. This seemingly contradictory result demonstrated that, prior to initiating therapy, PCC showed a small change in tumor size (less than 3E-5), so the negative correlation may not have a significant effect on the change in tumor size. PRCC showed this difference to be more exaggerated, indicating that it was essential to examine both PCC and PRCC. Changes in ν after 14 days of therapy were positively correlated with tumor size. This relationship was due to the negative feedback by which the increase in ν caused by T cells leads to an increase in PD-L1, and a decrease in T cells by increasing the amount of PD-1 and PD-L1 binding. This mechanism increased tumor size and suggested that high PD-L1 expression levels in T cells could inhibit the efficiency of IT. Likewise, an increase in expression level ν · µ expressed on tumor cells caused an increase in PD-L1 regardless of therapy. Additionally, an increase in ν · µ resulted in an increase in tumor size for similar reasons.

Discussion and Conclusions
Localized IR has been mediated tumor reduction in a T cell-dependent fashion. Therapeutic blockade of the T cell negative regulator PD-L1 could enhance T cell effector function when PD-L1 was expressed in tumors, as it happens in the tumor microenvironment after IR [22,23]. PD-L1 expression in the tumor microenvironment provided an opportunity for therapeutic intervention using immunotherapies such as anti PD-L1, anti PD-1 and CTLA-4 [24,25].
The purpose of this work was to obtain insight into mechanisms of the tumor-immune interactions while combination therapy of IR and IT was conducted. Notably, we considered PD-1 and PD-L1 in the model to discuss how these targets change. Mathematical modeling is one of the most important tools in analyzing the characteristics of the synergistic combinations and providing some useful insights about the dynamic of the interactions. A mathematical model for the tumor-immune interactions was developed. QSSA was used to simplify the model. Additionally, two-compartment models for treatment were considered to reflect the delay of tumor repression. The simulation was implemented for control, IT only, and IR only data. From the estimated parameters, combination therapy with IR and IT was fit. Simulations of our model are shown to be in reliable agreement with mouse experiments for all cases [14].
The model addressed how tumors and T cells were inhibited or stimulated by changes in PD-1 and PD-L1 expression. To this end, we utilized global sensitivity analysis to examine the influence of four parameters and confirmed the corresponding changes in tumor size. It was possible to employ a model to investigate how T cells and tumor growth were stimulated by PD-L1 and PD-1. The model predicted that T cells would have a significant influence on tumor growth and inhibition based on the ratio of PD-1 to PD-L1 when IT is applied. This result was explored in the global sensitivity analysis of the parameter changes, rather than with local changes. We compared PCC and PRCC from LHS among global sensitivities. PCC was intuitive, but the relationship between parameters and tumor size variation was more apparent when assessed with PRCC than PCC.
Our approach complemented Deng et al.'s experiments that combination therapy of IR and IT was likely an important part of the complex regulation in the IR tumor microenvironment that suppresses antitumor immunity. Two models can be contrasted. The model of Nikolopoulou et al. was a system of three ODEs with 19 parameters. They theoretically investigated the change in tumors against anti-PD-1. Their major contribution was to predict the general behavior of the tumor-immune response to therapy. Sensitivity analysis using PRCC via LHS and bifurcation analysis for tumor and T cells were also conducted. The model of Chappell et al. was a system of 4 ODEs with 16 parameters. Their study concerned tumor-T cell interactions with combination therapy of IR and IT. Regardless of the consideration of PD-1 and PD-L1, one of their contributions was reproducing the synergistic effects of IR and IT. Theoretical analysis such as model equilibria and stability was also conducted. Our approach reflected those models and additionally focus on the change in the expression levels of PD-1 and PD-L1. However, the proposed model has some limitations. We assumed that PD-L1 is permanently expressed on tumor cells and the expression increases with a tumor progression. However, PD-L1 is heterogeneously expressed on cancer cells and can be decreased over a period. Furthermore, we did not consider other immune checkpoint molecules upregulated on cancer cells that also have a significant impact on the cancer cell immune escape.
In Figures 6-9, the ratio rate η determined the ratio of PD-1 to PD-L1 in T cells. Before therapy, η and tumor size were positively correlated. That is, in our model, after cancer cells were implanted and before therapy, the binding of PD-1 to PD-L1 was small, and as a result, T cell suppression was small and the increase in tumor size was minimal, as shown in the experimental data ( Figure 6). The linearity was deregulated after beginning combination therapy (Figures 7-9). IT depleted PD-L1 to reduce the amount of binding to PD-1, thereby reducing T cell suppression and consequently suppressing tumors. IR additionally eliminated tumors and reduced the amount of PD-L1. Therefore, tumor size was determined based on T cell decrease/increase, where a decrease was caused by IR directly or an increase was caused by inhibiting the binding of PD-1 to PD-L1 indirectly.
Interestingly, as η decreased, the tumor size variation became more severe. The smaller the differentiation by T cells to PD-L1 was, the smaller the amount of PD-L1 was, and thus tumor growth was suppressed by IT. Meanwhile, when η was small with sufficient PD-L1, T cells increased the expression level of PD-1, thereby increasing the concentration of the PD-1-PD-L1 complex. As a result, the proliferation of the T cells was suppressed, and the tumor size increased. Therefore, two opposite cases can occur. Combination therapy overcomes the challenge by inhibiting the tumor both directly and indirectly. Using this model, we examined ways to regulate PD-L1 and PD-1 to improve the efficacy of combination therapy. After reducing PD-L1 and PD-1 expression by inhibiting tumor and T cells through IR, further depletion of PD-L1 through IT resulted in a decrease in the amount of PD-1-PD-L1 complex, thereby preventing the inhibition of T cell growth. This result indicated that if PD-L1 is lowly expressed in T cells, the efficacy of IT improves. That is, the level of PD-1-PD-L1 complex can be lower by reducing η.
It was noted that when η was small (i.e., when PD-1 was highly expressed in T cells or PD-L1 expression was low in T cells) and PD-L1 expression on the tumor was high, the tumor could be less inhibited as the PD-1-PD-L1 complex concentration was increased, indicating that T cell growth may be further suppressed. Additionally, combination therapy should be considered in two cases: when PD-L1 expression is high and when PD-1 expression is low or high in T cells. The expression level of PD-L1 was determined based on combination therapy. That is, when the expression level of PD-1 in T cells is low, successful therapy can be achieved by suppressing PD-L1 alone. Still, in the case of high expression of PD-L1 in T cells, such a strategy may fail due to excessive T cell reduction, although the tumor size was also reduced.
We investigated the tumor-immune interaction through a mathematical model. Through combination therapy, we verified the synergistic effect of combined IR and IT after parameter estimation through data fitting. Our model analyzed the mechanism of tumor inhibition through the relationship between PD-1 and PD-L1. We also examined the influence of each parameter through global sensitivity analysis of crucial parameters associated with the change in PD-1 and PD-L1. However, this analysis assumed a uniform distribution considering a slight variation in the parameters. Importantly, this assumption could not suggest a feasible condition for the parameters. In the next study, we will provide the distribution of the parameter space using Bayesian inference. This study could be valuable in that the Bayesian inference may be used for consideration of multiple doses and combination therapy remains challenging in the tumor-immune system.