Multi-Level Computational Modeling of Anti-Cancer Dendritic Cell Vaccination Utilized to Select Molecular Targets for Therapy Optimization

Dendritic cells (DCs) can be used for therapeutic vaccination against cancer. The success of this therapy depends on efficient tumor-antigen presentation to cytotoxic T lymphocytes (CTLs) and the induction of durable CTL responses by the DCs. Therefore, simulation of such a biological system by computational modeling is appealing because it can improve our understanding of the molecular mechanisms underlying CTL induction by DCs and help identify new strategies to improve therapeutic DC vaccination for cancer. Here, we developed a multi-level model accounting for the life cycle of DCs during anti-cancer immunotherapy. Specifically, the model is composed of three parts representing different stages of DC immunotherapy – the spreading and bio-distribution of intravenously injected DCs in human organs, the biochemical reactions regulating the DCs’ maturation and activation, and DC-mediated activation of CTLs. We calibrated the model using quantitative experimental data that account for the activation of key molecular circuits within DCs, the bio-distribution of DCs in the body, and the interaction between DCs and T cells. We showed how such a data-driven model can be exploited in combination with sensitivity analysis and model simulations to identify targets for enhancing anti-cancer DC vaccination. Since other previous works show how modeling improves therapy schedules and DC dosage, we here focused on the molecular optimization of the therapy. In line with this, we simulated the effect in DC vaccination of the concerted modulation of combined intracellular regulatory processes and proposed several possibilities that can enhance DC-mediated immunogenicity. Taken together, we present a comprehensive time-resolved multi-level model for studying DC vaccination in melanoma. Although the model is not intended for personalized patient therapy, it could be used as a tool for identifying molecular targets for optimizing DC-based therapy for cancer, which ultimately should be tested in in vitro and in vivo experiments.


Model construction
Based on the model scheme given in Figure 1 we developed a multi-level model accounting for different stages of DC-based anticancer immunotherapy. The model is composed of three parts representing different stages of the treatment -DC distribution in the human organs after injection, the biochemical pathway that regulates DC maturation, and the DC-mediated CD8 + T-cell responses.
Based on previous work (Ludewig, et al. 2004) we developed a mathematical model based on ordinary differential equations (ODEs) that accounts for the trafficking and distribution of DCs in the human body. Specifically, the injected cells mainly spread into the liver, the lung, the spleen, and other peripheries (Ludewig, et al. 2004, Mackensen, et al. 1999).
Furthermore, we describe the DC maturation by the activation of the biochemical signaling pathway NF-κB. The model is also based on ODEs and adopted from our previous work (Schulz, et al. 2017).
NF-κB is regulated by inhibitory proteins such as IκB. Stimulation with TNF-α, CD40L, and IL-1β leads to degradation of IκB what results in NF-κB activation. Since mature DCs secrete a variety of cytokines it was decided to concentrate on the most important ones in respect of T cell activation. IL-12 promotes IFN-γ production by T lymphocytes but it can also induce their proliferation (Henry, et al. 2008, Vecchio, et al. 2007). IFN-γ has antiviral, immunoregulatory, and anti-tumor functions (Parker, Rautela und Hertzog 2016, Tannenbaum und Hamilton 2000. The expansion and activation of T cells are also regulated by IL-6 which has a broad effect on immune cells (Hunter und Jones 2017). The chemokine IL-8 is included in the model because of its role in recruiting T helper cells (Taub, et al. 1996). DCs also express different surface molecules after stimulation. One of those is CD70. It supports the expansion of antigen-specific CD8 + T cells (Arens, et al. 2004, Borst, Hendriks und Xiao 2005. The NF-κB pathway activates through the stimulation of TNFR1 and IL-1R receptors by TNF-α and IL-1β, respectively. CD40L activates the pathway through the CD40 receptor. After recognition of those proteins, recruitment of protein effectors can result in phosphorylation. For the model validation, it was necessary to add a stimulation with LPS which activates the TLR-4 receptor. The last step of our model accounts for the DC-mediated activation of CD8 + T cells in the spleen. An adult human contains approximately around (3-5)⋅10 11 naive T cells (Abbas, Lichtman und Pillai 2018, Jenkins, et al. 2010, whereof one cell out of 10 5 -10 6 expresses a sufficient TCR for a specific antigen (Celli, et al. 2012, Henrickson, et al. 2008. After the naive T cells have matured in the thymus they continuously move through the blood, lymphatic vessels, secondary lymphoid organs like the spleen and non-lymphoid tissue for antigen recognition. After homing into the spleen, they spend around 24 h in the T cell zones, where the interaction with DCs occurs (Celli, et al. 2012, Nolz 2015. The T cell activation happens in three phases (Henrickson, et al. 2008, Ozga, et al. 2016: the first phase is a short-term interaction between the naive T cell and the DC, it lasts about 8h and the T cell starts to upregulate activation markers. The second phase, which lasts about 12 h, is accompanied by further upregulation of activation markers and initiation of IFNγ and IL-2. The third phase begins about 24 h after the first contact, during this short T cell-DC interaction the T cell proliferation is induced. In addition, it was assumed that the activated T cells also remain in the spleen for about 24 h, as they usually migrate to the site of infection immediately after activation. In addition, the mortality rate of the naive T cells can be neglected since they have a half-life of several years (Jenkins, et al. 2010, Vrisekoop, et al. 2008 and in this study, only a period of a few weeks is considered. It is not yet fully understood whether T cells follow a linear or asymmetric differentiation pathway (Ahmed, et al. 2009). Nevertheless, this work will follow an asymmetric approach as proposed in several other works and models (Arsenio, et al. 2013, Crauste, et al. 2017, Joshi, et al. 2007, Girel, et al. 2007). Thus, DC-induced naive T cell activation leads to differentiation into early effector CD8 T cells. Those cells proliferate after the second contact with antigen-presenting cells as described in the above paragraph. After some time, the early effectors differentiate into short-lived effector cells and memory precursor cells generating long-lived memory cells. The process of whether the early effector cells become short-lived effector or memory precursor cells is supported by the available amount of inflammatory cytokines like IL-12 during priming, i.e. high inflammation induces the generation of more short-lived effector cells (Joshi, et al. 2007, Pearce und Shen 2007).
Short-lived effector CD8 T cells degrade with a half-life of approximately 65-80 days (Joshi, et al. 2007), thereby around 10% become long-lived memory cells (Mueller, et al. 2013, Badovinac, Porter und Harty 2002. In addition, two different types of long-lived memory T cells can be distinguished: effector memory T cells and central memory T cells (Arsenio, et al. 2013, Wherry, et al. 2003. It is still controversial under which conditions effector memory T cells and central memory T cells occur (Ahmed, et al. 2009), thus there is no reliable data available. Therefore, memory precursors and all kinds of long-lived memory cells are generally referred to as memory T cells within this work.
In the following, the model equations are presented. For better readability, the model is divided into the three main components in the paragraphs below.

DC distribution
The following system of ordinary differential equations (ODE) describes the distribution of DCs in the human body after injection. We describe the temporal cell dynamics in the blood, the spleen, the liver, and the lung.
Equation (1) accounts for gradual decrease of the DCs in the blood ( ).
Where the parameter describes the total distribution rate into the other compartments. denotes the homing rate into other periphery.
Equation (2) shows the temporal evolution of DCs in the spleen ( ).
Thereby the cells have a homing rate of from blood into the spleen and leave the compartment with a rate of 0 . The splenic volume is given by .
The DC distribution into the liver ( ) is given by equation (3). The amount of cells remains constant over time (Mackensen, et al. 1999), while the homing rate from blood to liver is given by . The volume of the liver is denoted by Equation (4) displays the DC activity in the lung ( ). It was found that the cell number does not shrink below a certain threshold denoted by (Mackensen, et al. 1999). The homing rate is given by whereby the lung is cleared again with a rate of 0 and is the volume of the lung. The initial concentrations of DC in the spleen, lung and liver were assumed to be zero.
All four ordinary differential equations can be solved analytically by separation of variables or variation of constants if the equation is nonlinear. Thereby the following system of equations is obtained: With the initial conditions chosen as follows:

DC activation
After stimulation with TNF-α, IL-1β, CD40L, or LPS the respective receptor-associated proteins are phosphorylated which leads to activation of IKK. Thereby IκBα is degraded which results in the release of free NF-κB. The transcription factor NF-кB enters the nucleus where the transcription of many genes is upregulated. Among these genes, the model accounts for IκBα and the most important ones for T-cell activation. This further leads to the secretion of chemokines and cytokines, and the upregulation of surface molecules, i.e., IL-12p70 mRNA, IL-6 mRNA, IL-8 mRNA, and CD70 mRNA.
The phosphorylation of IRAK1 and TRAF2 is described as follows: Stimulation with IL1-β and LPS lead to phosphorylation of IRAK1 with the rates ℎ1 1 and ℎ2 1 , respectively. The synthesis and degradation rate of IRAK1 are denoted by 1 and 1 , respectively. IRAK1p describes the phosphorylated state with a degradation rate of 1 .
This is equivalent for the phosphorylation of TRAF2 by the stimulation with TNF-α and CD40L. Here the phosphorylated state is denoted with TRAF2p.
This process leads to a state transition of IKK: Where is the synthesis rate of IKK, 1 and 2 are the respective activation rates of IKK through the phosphorylated proteins IRAK1p and TRAF2p. The degradation rate is given by .
Activation of IKK results in IKKβ increase, which is degraded with a rate : This activation can release free NF-κB by degrading IκBα: is the gain rate of free NF-κB in the nucleus mediated by IκBα and IKKβ and is the loss of free NF-κB mediated by IκBα. The total amount of free NF-κB is given by Ntot. 1 and 2 are the Michaelis-Menten like coefficients.
The transcription factor NF-κB enters the nucleus and upregulates the transcription of many genes, including the transcription of gene for IκBα into mRNA (mIκBα): With transcription rate and degradation rate .
This mRNA is further translated into IκBα which again leads to down regulation of NF-κB: denotes the translation of mIκBα into IκBα whereby the is the loss of free IκBα mediated by NF-κB and IKKβ. 3 is the Michaelis-Menten like coefficient.
The transcription of mIL-8 and mIL-6 is also upregulated: are the respective transcription rates induced by NF-κB and the degradation rates of mIL-8 and mIL-6 are given by 8 and 6 .
Those are translated into the chemokine IL-8 and the pro-inflammatory cytokines IL-6, respectively: With the translation rates 8 and 6 and the degradation rates 8 and 6 . The secretion of the proteins by the DCs is mediated by a rate of 8 for IL-8 and 6 for IL-6.
NF-κB also stimulates mCD70 which results in a presentation of CD70 molecule on the cell surface: Again 1 70 describes a basal transcription rate, whereby 2 70 is the rate of transcription of mCD70 mediated by NF-κB.
70 is the translation of mCD70 into CD70 and 70 is the degradation of the surface marker.
The cytokine IL-12p70 is not directly activated because it is a heterodimer consisting of IL-12p40 and IL-12p35 joined by bond. IL-12p40 is constantly active but IL-12p35 is activated through NF-κB.
For simplicity it is estimated that the amount of IL-12p40 is very small compared to IL-12p35. So its activation can be written as a direct result in IL-12p70 translation: and degradation rate 12 35 of mIL-12p35. As well as the translation rate of mIL-12p35 into IL-12p70 12 35 , the degredation rate 12 70 and the secretion rate 12 70 of IL-12p70.
The secretion of the cytokines can be described as follows: Where 12 70 , 6 and 8 are the distribution rates for IL-12p70, IL-6, and IL-8 in blood, respectively.
The initial amount of the variables was chosen to be zero, except for IRAK1(t = 0), TRAF2(t = 0), IKK(t = 0), and IκBα(t = 0). Those initial values were set equal to one. If activation through Lipopolysaccharide (LPS) is necessary the respective parameter has to be set equal to one, whereof the other input signals TNFα, CD40L, and IL1-β have to be zero. For activation through TNFα, CD40L, and IL1-β these values have to be chosen as one and LPS equal to zero.

DC-mediated T-cell response
The last part of the model describes the dynamics of T cell activation through contact with DCs.
Thereby we concentrate on the spleen, as we currently have no data available for other lymphoid organs. Entering the spleen DCs travel to T-cell zones where the interaction takes place. After successful contact between DCs and naive T cells, the T cells are activated and differentiate into early effector T cells. After second contact, early effector cells proliferate and further differentiate into short-lived effector cells and long-lived memory cells. Whereof short-lived effector cells have a half-life of approximately 65-80 days and around 10 % become memory T cells (Mueller, et al. 2013, Badovinac, Porter und Harty 2002. We will neglect the homing of activated short-lived effector and memory T cells into other compartments or the tumor microenvironment since we are mainly interested in the total amount of memory cells.
All together results in the following system of differential equations accounting for the DC-T cell interaction: Equation (33) describes the activation of naive T cells ( ). Thereby, the cells home into the spleen with a rate of ℎ . 0 is the initial value for naive cells that get activated with a rate of . The activation depends on the amount of available DCs in the spleen ( ) and the intensity of cytokines, chemokines, and surface markers produced or presented by DCs characterized by the following function Equation (37)  The temporal dynamics of is characterized by equation (35), where is the degradation rate.
A fraction (10%) of degrading becomes memory T cells.
Equation ( The initial conditions of the model variables are given as follows:

Parameter estimation for different parts of the model
To calibrate the model, we used experimental data to adjust the parameter values. The optimization strategy is described in the Materials and Methods. For better readability, this section is divided into three parts.

DC distribution
For calibration, the model simulation was fitted to experimental data from Mackensen et al. (Mackensen, et al. 1999). The data describes the distribution of in vitro differentiated DCs into the liver, the spleen, the lung and other periphery after intravenous injection. The result is given in Figure 2.

DC activation
Since it was not possible to find enough data to fit the model some parameters were kept constant to reduce their number and to prevent overfitting. Synthesis rates were set equal to the degradation rates to make the steady states equal to the initial values when no input to the system is given. Further, the translation constants 1 12 70 , 1 6 and 1 8 could be set to zero because the proteins are not basal as shown by Pfeiffer et al. (Pfeiffer, et al. 2014). The values for 1 , 2 and 3 were set equal to one for simplicity. Since so far, no dynamical data for the NF-κB pathway activation in dendritic cells with activation through TNF-α, CD40L, or IL-1β was found, but with LPS. Therefore, we used the experimental data from Bode et al. (Bode, Schmitz, et al. 2008) to calibrate the model and used LPS as only the input signal. The data contain information about the dynamics of IκBα mRNA, IκBα, and NF-κB. For the IκBα fitting, additional simulation points were added since the given data did not account for the entire period being simulated. The result for the model simulation compared to the data is given in Figure 3. To validate the dynamics for the cytokines IL-12p10 and IL-6, the chemokine IL-8, and the surface marker CD70 we used the data set from Pfeiffer et al. (Pfeiffer, et al. 2014). Therefore, TNF-α, CD40L, and IL-1β but no LPS were used as input signals. The outcome is shown in Figure 4. The concentrations were measured 24 h after electroporation with caIKKβ-RNA.
To simulate electroporated DCs we introduced a modified degradation rate of IKKβ, i.e., , .
We also adjusted this parameter value to calibrate the model.

DC-mediated T-cell responses
To calibrate the DC-T cell interaction part of the model experimental data from Pfeiffer et al. (Pfeiffer, et al. 2014) is used. In this work 1 ⋅ 10 6 CD8 + T cells are primed in vitro with Mock-electroporated and caIKKα-and caIKKβ-mRNA (alone or in combination) electroporated DCs. The DCs were used four hours after electroporation to stimulate the T cells at a ratio 1:10. After one week the total number of MelA-specific CD8 + T cells was determined. The T cells were restimulated two times. Those caIKKα-and caIKKβ-mRNA electroporated DCs have up to 10 times higher cytokine/chemokine production and can stimulate nearly five times more T cells (Pfeiffer, et al. 2014).
For the model calibration caIKKβ-RNA electroporated DCs are simulated by the modified degradation rate of IKKβ, i.e., , . To simulate the experimental conditions in the model, the aboveintroduced system of equations is reduced to: With the number of DCs fixed to = 1 ⋅ 10 5 . To fit the reduced model to the data, the parameter 4 was optimized. The results of the model calibration is given in Figure 5. Both model variations were optimized in parallel.

Sensitivity Analysis
We performed global sensitivity analyses to quantify the impact of model parameters on the dynamics of the system. We used the Sobol method that considers variations within the entire variability space of the model parameters (Saltelli, et al. 2008, Sarrazin, Pianosi und Wagener 2016. The analysis was performed using the MATLAB toolbox SAFE (Pianosi, Beven, et al. 2016).
Specifically, the contribution to the variance of model output (e.g., a variable accounting for the number of immune cells) by an individual input factor (i.e., a model parameter denoting the degradation of a protein) was considered a good measure for sensitivity in our analysis. Considering the following definition for the model response with the model output, ⃗ = ( 1 , … ) the vector containing the input factors, the number of input factors, and a function mapping input to output, there are different ways to calculate this contribution: first-order indices (main-effects) and total-order indices (total-effects) are considered.
By definition (Saltelli, et al. 2008, Sarrazin, Pianosi und Wagener 2016, the first-order indices are calculated using the following equation and the total-indices are defined as follows Where is the i-th input factor and ∼ is the vector of all input factors except the i-th one. We generated the input samples using the Latin-Hypercube sampling method that not only samples random parameter values but also guarantees a uniform distribution of parameter values in their defined ranges (Sarrazin, Pianosi und Wagener 2016). Assuming a random sampling of size thereby leads to = ( + 2) model evaluations due to the resampling strategy (Sarrazin, Pianosi und Wagener 2016). We chose to be 2500.
To assess the uncertainty of the global sensitivity analysis, we investigated the convergence of the sensitivity indices using the SAFE toolbox. A value of the width of the confidence interval (here 95%) close to zero indicates that the sensitivity index converges (Pianosi, Beven, et al. 2016). The maximum width of the confidence interval was calculated by where and are the upper and lower bounds of the sensitivity index of the i-th input factor.     (Mackensen, et al. 1999) Human experiment volume liver 0.2 ⋅ fixed - (Sauermost und Freudig 1999) Human data volume spleen 0.16 ⋅ fixed - (Sauermost und Freudig 1999) Human data volume lung 0.3 ⋅ fixed - (Sauermost und Freudig 1999) Human data volume blood 6-7 liter -- (Sauermost und Freudig 1999) Human data Number of DCs injected in blood 1 ⋅ 10 7 cells fixed - (Mackensen, et al. 1999) Human data number of DCs that stay in lung 0.13 ⋅ 10 7 cells 0.13 ⋅ - (Mackensen, et al. 1999) Human data