We developed a multi-level model accounting for three stages of DC immunotherapy – the spreading and distribution of intravenously injected DCs in the human organs, the biochemical reactions regulating the DCs’ maturation and activation, and DC-mediated T-cell responses. The model was implemented using ODEs and simulated in MATLAB R2015b (see Supplementary Material for details). We used the MATLAB function ode45 to solve the system accounting for the maturation of DCs. The function is based on the Dormand and Prince Runge-Kutta methods. We used the MATLAB function dde23 to solve the equations accounting for a DC-mediated T-cell response. This method is based on an explicit Runge-Kutta method pair and especially for delay differential equations with constant delay. We ran simulations on a computer with a four core CPU (3.2 GHz) and 8 GB RAM.

Here, we listed some representative equations for each part of the model. Specifically, we used four ODEs to describe the temporal dynamics of DCs in the blood, the spleen, the liver, and the lung. A representative equation is shown as follows d d t D C S p l e e n ( t ) = μ B S Q B l o o d Q S p l e e n D C B l o o d ( t ) − μ S 0 D C S p l e e n ( t ) (1)

The number of DCs in the spleen ( D C S p l e e n ) is determined by the flow of DCs from the blood ( D C B l o o d ) into the spleen at the rate ( μ B S ) and its degradation rate ( μ S 0 ). Q B l o o d and Q S p l e e n denote the volume of the blood and spleen, respectively.

The maturation of DCs in the spleen is characterized by 20 ODEs that account for the activation of the NF-κB pathway and its downstream targets such as cytokines that are crucial for T-cell responses. The equations accounting for the mRNA and protein of IL-8 are shown below d d t m I L 8 ( t ) = k t r a n s c 1 m I L 8 + k t r a n c s 2 m I L 8 ⋅ N F κ B ( t ) − k d e g m I L 8 ⋅ m I L 8 ( t ) (2) d d t I L 8 D C ( t ) = k t r a n s l I L 8 ⋅ m I L 8 ( t ) − k d e g I L 8 ⋅ I L 8 D C ( t ) − k s e c I L 8 ⋅ I L 8 D C ( t ) (3)

The transcription of IL-8 mRNA is determined by its basal transcription rate ( k t r a n s c 1 m I L 8 ) and another term ( k t r a n c s 2 m I L 8 ⋅ N F κ B ( t ) ) accounting for the regulation by NF-κB. k d e g m I L 8 denotes the degradation of the IL-8 mRNA. The protein expression of IL-8 is determined by its translation from its mRNA ( k t r a n s l I L 8 ), its degradation ( k d e g I L 8 ), and its secretion ( k s e c I L 8 ) from DCs into the spleen.

The activation of T cells is modeled using four ODEs. The equations describe the process of how lymph node T cell activation translates to different phenotypes of T cell subsets such as early effector, short-lived effector, and memory T cells. The equation describing memory T cells is shown as follows d d t M ( t ) = k d i f f 2 E E ⋅ E E ( t ) + 0.1 ⋅ k d e g S L E ⋅ S L E ( t ) (4)

The number of memory T cells (M) is translated from early effector (EE) and short-lived effector (SLE) T cells with the rates of k d i f f 2 E E and k d e g S L E , respectively. The constant 0.1 denotes that in experiments only 10% of the short-lived effector T cells become memory T cells (Badovinac et al., 2002; Mueller et al., 2013).

Before model calibration, it is useful to investigate whether it is possible to obtain identifiable parameters using experimental data. Global structural identifiability analysis can provide a good indication of this. Therefore, we used the MATLAB toolbox GenSSI to perform structural identifiability analysis (Chis et al., 2011). The algorithm uses Lie derivatives of the ODE model to investigate the structural identifiability of ODE models. Theoretically, if sub-models are structurally identifiable, so is the entire model (Villaverde et al., 2016). Therefore, to reduce computation time we divided our model into separated parts, including DC distribution, NF-κB activation, NF-κB-mediated secretion of cytokines, and T-cell response, and performed the analysis on each part.

We used a hybrid method that combines global and local optimization algorithms to perform parameter optimization. Such a method facilitates global exploration of parameter space and fast local convergence (Villaverde et al., 2019). Specifically, we derived 1000 parameter sets using Latin hypercube sampling that not only samples random parameter values but also guarantees a uniformed distribution of parameter values in their defined ranges (Tang, 1993). The 1000 parameter sets are initial values of model parameters and were used for model calibration. We first fit the model to the experimental data using the pattern search algorithm (MATLAB function patternsearch) that is a global derivative-free optimization algorithm. The top 100 the solutions of the global optimization results (quantified by the cost function) were used for subsequent local optimization. The local optimization algorithm (MATLAB function fmincon) is gradient-based and allows for efficient searching that makes the cost function converge fast. We obtained the optimum parameter set that minimizes the cost function Φ ( p ) = ∑ i , j ( y i ( t j ) − y i ( t j , p ) m a x j ( y i ( t j ) ) ⋅ s d ( y i ( t j ) ) ) 2 , (5) Where y i ( t j ) represents the experimental data that shows the value of the observation y i ( t j ) at time point j. y i ( t j , p ) is the corresponding model simulation with a specific set of parameter values p. The cost function is normalized using the maximum value of each experimental data set (i.e., m a x j ( y i ( t j ) ) ) to prevent the biased effects caused by different data scaling during parameter estimation. If available, the cost function is additionally weighted by the standard deviation of the experimental data s d ( y i ( t j ) ) .

We used the experimental data accounting for in vitro differentiated DCs into the liver, the spleen, the lung and other periphery after intravenous injection (Mackensen et al., 1999) to characterize the dynamics of DCs in human organs (Supplementary Table S1). We characterized NF-κB pathway activation in DCs (Supplementary Table S2), cytokine and chemokine production by DCs (Supplementary Table S2), and DC-mediated T cell responses (Supplementary Table S3) using the data from our previous publication (Pfeiffer et al., 2014). A detailed description of how each part of the model is calibrated can be found in Supplementary Material.

To analyze the uncertainty in parameter estimates, we computed Pearson correlation coefficients to quantify their linear dependence using the top 100 out of 1000 estimates. Among the top 100 parameter estimates, we used the best 15 parameter estimates that show the minimum value in the cost function to calculate the confidence intervals of estimated parameters. Specifically, for each estimated parameter we generated 1000 bootstrap samples using its estimated values in the best 15 parameter estimates. Then, we used the mean value and standard deviation of the 1000 samples to derive the parameters’ 95% confidence intervals. We performed the analyses using the MATLAB function bootci.

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 et al., 2016). We computed two types of indices: first-order indices (main-effects) and total-order indices (total-effects). The former measures the direct contribution from a model parameter (e.g., the total amount of NF-κB) to a model variable (e.g., the count of memory T cells), while the latter measures the overall contribution including the direct contribution and the amplification of this contribution due to interactions with all other model parameters (Sarrazin et al., 2016). The analysis was performed using the MATLAB toolbox SAFE (Pianosi et al., 2016) and the detailed computation of the sensitivity indices can be found in Supplementary Material.

We developed a multi-scale model accounting for DC-based anti-cancer immunotherapy. We considered different stages of the DC therapy: 1) the bio-distribution of DCs in the human body after the treatment, 2) the biochemical pathways underlying DC maturation, and 3) the DC-induced immune response such as activation of CD8⁺ T cells.

We characterized the trafficking and distribution of dendritic cells (DCs) in the human body using published data (Ludewig et al., 2004) (Figure 1; see Supplementary Material for details). Specifically, DCs are administrated into patients through intravenous injection. The injected DCs mainly spread into the liver, the lung, the spleen, and other peripheries. According to the data (Mackensen et al., 1999; Ludewig et al., 2004), after reaching the liver, DCs reside there, while DCs enter quickly into the lung but also decrease to a minimal level. In our model, we used the spleen as a representative lymphoid organ, in which effective, antigen presentation-mediated interactions between T and dendritic cells happen (Mackensen et al., 1999; Barinov et al., 2017; Abbas et al., 2018).

Concerning the pathways underlying DC maturation and activation, we further considered in the model the NF-κB pathway underlying DC maturation; this part of the model was adapted from our previous results (Schulz et al., 2017) (Figure 1; see Supplementary Material for details). In our model, mature DCs secrete a variety of cytokines (e.g., IL-12, IL-6, and IL-8) that can lead to T cell activation. The production of the cytokines is due to the stimulation of receptors (such as TNFα receptor, IL-1 receptor, CD40 receptor, and TLR4), leading to the activation of the NF-κB signaling pathway. Besides, we considered a surface protein CD70 expressed by DCs after stimulation, as the protein can induce the expansion of antigen-specific CD8⁺ T cells.

Finally, we included a model module accounting for T cell activation by DCs in the spleen. The process was modeled considering three phases: 1) a short-term interaction between the naive T cells and DCs; 2) upregulation of activation markers and initiation of IFNγ and IL-2; and 3) T cell proliferation after contact with DCs.

Taken together, we developed a model accounting for the life cycle of intravenously injected DCs from their spreading in human organs to inducing a T-cell response in the spleen. The resulting model contains 25 variables and 46 parameters (see Supplementary Material for details). The model includes not only cell population dynamics in human organs but also biochemical reactions that are crucial for DC maturation, making it an in silico platform to investigate intracellular manipulation of DCs utilized in vaccination against the tumor.

After constructing the multi-scale model, we performed structural identifiability analysis to identify whether model parametrizations with different values can produce the same simulations results (see Materials and Methods). The results showed that the whole model is not structurally identifiable. Specifically, the DC distribution part of the model is structurally non-identifiable, and this is caused by the parameters accounting for volumes of human organs (Q _Blood , Q _Spleen , Q _Lung , and Q _Liver ), So, we fixed their values using the physiological information (Supplementary Table S1). The part accounting for NF-κB activation is structurally non-identifiable, and this is caused by parameters accounting for the phosphorylation rate of IRAK1 ( k p h 2 I R A K 1 ) and TRAF2 ( k p h 2 T R A F 2 ). Hence, we made their values equal to other phosphorylation processes catalyzed by other enzymes (Supplementary Table S2). The other model parts are structurally identifiable, and they account for NF-κB-mediated secretion of cytokines and chemokines and T cell response.

Next, we characterized the model parameters using published experimental data sets (see Materials and Methods). We separately calibrated the model using independent datasets that account for the dynamics of the system at different levels. Such a strategy is suitable and has been used for biological models composed of parts with different structural, time, and space scales (Vera et al., 2013).

We used the data that quantify injected DCs’ activities in the lung, liver, spleen, and blood to characterize model parameters associated with the trafficking and distribution of DCs in the human body (Mackensen et al., 1999). The obtained model can reproduce the data available (Figure 2). After DCs are injected into the blood, they quickly spread into the other organs and the amount of DCs decreases to zero in the blood. DCs traffic into the lung leading to a temporal increase followed by a quick drop, and afterward, the DCs stay at a low level. DCs enter the liver and remain stable in number for up to 72 h (Mackensen et al., 1999). In the spleen, the amount of DCs rapidly reaches a peak and gradually decreases as they circulate in the body after provoking a T-cell response.

As the model equations accounting for the NF-κB pathway underlying DC maturation were adapted from our previous model, we used their original parameter values as initial values and re-calibrated them by fitting the model simulations to the experimental data measuring NF-κB pathway activation after LPS stimulation (Bode et al., 2009). Using the data, we characterized the dynamics of IκBα mRNA and protein and NF-κB in DCs (Figure 3). The LPS stimulation upregulates the expression of free NF-κB, as the increased IKKβ by the stimulation releases NF-κB from the complex formed by NF-κB and IκBα and degrades IκBα through phosphorylation. The free NF-κB promotes the transcription of the IκBα mRNA, leading to the recovery of IκBα that downregulates the expression of free NF-κB through a negative feedback loop. Besides, we used another set of data to characterize the dynamics of cytokines and surface markers after electroporating DCs with RNAs encoding constitutively active IKKα and -β (caIKK) that can activate the NF-κB pathway (Pfeiffer et al., 2014). Such treatment results in upregulation of IL-12, IL-6, IL-8, and CD70 for 72 h (Figure 4) and these markers are crucial for priming a T-cell response.

To characterize the model parameters that are associated with T-cell priming by DCs in the spleen, we used the in vitro data that show dynamics of the T-cell population after co-culturing them with control DCs or caIKK-DCs (Pfeiffer et al., 2014). Compared to the control DCs, the caIKK-DCs secret more cytokines such as IL-12, IL-8, IL-6, and TNF (Pfeiffer et al., 2014) and increase the production of T cells (Figure 5).

Furthermore, we performed practical identifiability analysis in parameter estimates (see Materials and Methods). This allowed us to examine whether the estimated parameters are practically identifiable – the estimated model parameters have unique values that fit model simulations to experimental data used for model calibration. As shown in Figure 6A, the estimated parameters for DC distribution have no correlation with each other suggesting the corresponding parameters are practically identifiable. This is confirmed by the distribution of estimated parameter values in the best 15 parameter estimates that show the minimum cost function value (Supplementary Figure S1). All estimated parameters for DC distribution have unique values for the best parameter estimates. In addition, the model parameters accounting for the NF-κB pathway underlying DC maturation show moderate correlations (Figures 6B,C). Among the 28 estimated parameters, five are practically non-identifiable and they are k p h 1 T R A F 2 , k d e g T R A F 2 p , k d e g m I L 6 , k t r a n s c 2 m I L 8 , and k d e g I L 8 (Supplementary Figure S2). The non-identifiable parameters show small variances in their estimated values (Supplementary Table S2). This is due to the relatively small number of experimental data that are available for parameter estimation (Raue et al., 2009). The Michaelis-Menten coefficient K ₄ is the only parameter estimated to fit model simulations to the data accounting for T-cell dynamics (Supplementary Table S3). The estimated value of K ₄ is practically identifiable for its unique value in the best parameter estimates (Supplementary Figure S3). Taken together, most model parameters are practically identifiable because of their unique estimated values in the best parameter estimates. The practically non-identifiable parameters have confidence intervals with small ranges, suggesting minor influences on their biological interpretability.

Due to the lack of a complete data set that can characterize the dynamics of DCs used as an anti-cancer therapy, we calibrated model modules separately using relevant datasets from different publications. In all cases, the datasets were produced using human material or relevant experimental models that are generally accepted in the context of DC vaccine development (Brossart et al., 2001). We think they are complementary because they characterize the dynamics of the DC vaccination at different levels. This strategy has been used in other data-driven computational models (Sobotta et al., 2017, Hesse et al., 2021, Fey et al., 2015).

After finishing calibrating the model with experimental data, we used it to simulate a DC-mediated T-cell response. Before running simulations, we set DC _in = 10⁵ as the data showed that about 10⁵ DCs are required for a T-cell response with a 70% probability (Celli et al., 2012). We set Q _Spleen = 10⁵ as the measured volume of a spleen is 10⁵ mm³ (Odorico et al., 1999). Besides, we assumed that the initial number of antigen-specific naive T cells in the spleen is 10⁶, so we set T ₀ = 10⁶ (Celli et al., 2012, Henrickson et al., 2008). We simulated the dynamics of different T cells in the spleen after injection with two different DC vaccines: normal DCs and DCs electroporated with caIKKβ-RNA (caIKK-DCs). In the simulations, both types of DCs were injected at t = 0 h with different degradation rates of IKKβ that are caused by caIKK in DCs ( k d e g I K K b = 0.840 h⁻¹ for normal DCs and k d e g , m o d I K K b = 0.216 h⁻¹ for caIKK DCs).

As shown in Figures 7A–D, the DC vaccines result in the loss of naive T cells that differentiate into early effector T cells, which show a quick increase after the DC stimulation. The early effector cells gradually differentiate into short-lived effector T cells and memory T cells, both of which saturate at high levels. Compared to the normal DC vaccine, the caIKK-DCs increase the levels of short-lived effector and memory T cells by about 7-fold, demonstrating the enhanced immunogenic potency of the caIKK-DCs. At the molecular level, such improved immunogenic potency is caused by the upregulated activation of NF-κB pathway by IKKβ in DCs (Figures 7E,F).

Successful activation of naive T cells depends mainly on successful and strong interaction with antigen-presenting DCs (Timmerman and Levy, 1999, Abbas et al., 2018). The goal of a DC vaccine is to increase the number of the resulting memory T cells that contribute to a rapid immune response upon reactivation and form a long-lasting immunity (Akondy et al., 2017, Ando et al., 2019, Ahmed und Gray, 1996). Therefore, we performed sensitivity analyses to investigate the molecular mechanisms that are crucial for regulating the differentiation from naive T cells into early effectors and thus into memory T cells. Specifically, we used the global sensitivity method Sobol to compute sensitivities of model parameters to the population of memory T cells over the simulation time interval [0, 200] h (see Material and Methods). As shown in Figures 8A,B, the top-ranking 15 parameters show similar patterns in their sensitivity indices – the values are low shortly after the DC stimulation, gradually increase to a higher level, and stays at the high level until the end of the simulations. Among the top-ranking parameters, the most influential ones on the production of memory T cells are k d e g I K K b , k a c t N , and N _tot that account for the degradation rate of IKKβ, the activation rate of naive T cells, and the total amount of free NF-κB. The less influential parameters are those associated with the degradation rate of IκBα mRNA ( k m R N A m I k B a ), the degradation rate of the IKK protein ( k d e g I K K , short-lived T-cell differentiation ( k d i f f 2 E E ), loss of free IκBα ( k l o s s I k B a ), and the production of IκBα mRNA and protein ( k t r a n s c m I k B a and k t r a n l I k B a ). The least influential parameters are the degradation rate and NF-κB-mediated transcription rate of IL-8 mRNA ( k d e g m I L 8 and k t r a n s c 2 m I L 8 ), the degradation rate and NF-κB-mediated transcription rate of IL-6 mRNA ( k d e g m I L 6 and k t r a n s c 2 m I L 6 ), the injected number of DCs (DC _in ), and the homing rate of DCs into spleen ( μ B S ). We obtained similar results while computing the sensitivities of model parameters to the total amount of memory T cells (computed by taking the integral over the simulation interval [0, 200] h) (Figure 8C). The only difference is that the DC emigration rate from blood to other organs (µ) replaces k d i f f 2 E E and appears as the least influencing parameters.

Further analysis showed that the length of simulation time and the variation of estimated parameter values have limited effects on model parameters’ sensitivity indices for the total amount of memory T cells over the simulation interval. Specifically, after extending the simulation time to 500 h, the most influential parameters (i.e., k d e g I K K b , N t o t , k a c t N , k d e g m I k B a , and k d e g I K K ) remain unchanged (Supplementary Table S5). In longer stimulation time, several parameters (i.e., µ, µ _BS , DC _in , k t r a n l I k B a , k d e g m I L 6 and k t r a n s c 2 m I L 6 ) become less influential and two parameters k d e g I L 6 (the degradation rate of IL-6) and K ₄ become more influential. After increasing parameter variations to 90% of their estimated values, most of the top 15 parameters remain in the list but have different ranking (Supplementary Table S6). The exceptions are µ _BS , k t r a n s c 2 m I L 6 , and k d e g m I L 8 that become less influential and drop out of the top 15 parameters. Besides, k d e g T R A F 2 p (degradation rate of TRAF2), k d e g I L 8 , and k d i f f 2 E E (differentiation rate of early effector T cells into memory T cells) become more influential and are new top 15 parameters.

Taken together, the results demonstrated the ability of the multi-scale model to differentiate the ability of different DC vaccines to stimulate a T-cell response and to reveal the molecular mechanisms that are important for CD8⁺ T-cell activation through sensitivity analysis.

After identifying influential parameters that can modulate the production of memory T cells, we ran simulations to predict how corresponding biological processes can change the dynamics of the memory T-cell population. From the 15 most influential parameters, we have selected five corresponding to specific genes that could be experimentally manipulated. Specifically, we perturbed those parameters in an interval that decreases and increases their estimated values by 10 folds (i.e., the estimated value × [0.1, 10]) and computed the steady state of memory T cells. As shown in Figure 9, the total amount of free NF-κB (N _tot ) and the degradation rate of IκBα mRNA ( k d e g m I k B a ) can positively affect the population of memory T cells. Decreasing the value of N _tot by 90% eliminates the T memory cells while increasing its value leads to an increased cell population. The cell population peaks when the value of N _tot increases by about 3–4 folds and slightly decreases when N _tot is at its maximum level. Such a phenomenon could be explained by the negative feedback loop formed by NF-κB and IκBα. Free NF-κB activates the transcription of IκBα, whose encoding protein reduces free NF-κB by forming complexes. Thus, when the level of IκBα protein reaches a certain threshold, free NF-κB starts decreasing leading to a reduced number of T cells. In contrast, the decreasing of k d e g m I k B a by 90% results in a slight reduction of T memory cells, and a 10-fold upregulation in the parameter value leads to about a 3-fold increase in the cell population. Biologically, the IκBα protein is an inhibitor of NF-κB and traps free NF-κB through forming complexes (Hayden and Matthew, 2008), so decreasing the IκBα mRNA results in the reduced level of the protein, thereby releasing more free NF-κB in DCs that is required for T-cell activation.

On the other hand, perturbation of the degradation rates of IKKβ ( k d e g I K K b ), IL-6 mRNA ( k d e g m I L 6 ), and IL-8 mRNA ( k d e g m I L 8 ) negatively affect the memory T-cell population. The dynamics of the cell population show similar patterns when the three parameters are perturbed in the specified interval – memory T cells are at the maximum level when the values of the parameters are reduced by 90% and at the minimum level when the values of the parameters increase by 10-folds. Biologically, IKKβ induces degradation of IκBα through phosphorylation (Hayden and Matthew, 2008), thereby releasing NF-κB from the complexes. Therefore, increasing the degradation of IKKβ leads to downregulation of NF-κB in DCs that reduce the induction of memory T cells. IL-6 and IL-8 secreted by DCs are required for differentiation of naive T cells into early effector T cells that further differentiate into short-lived effector T cells and memory T cells (Hunter and Jones, 2017, Taub et al., 1996), thereby increasing the degradation of the cytokines results in the decreased level of memory T cells.

Furthermore, we simulated how the population of memory T cells changes when combining two parameters and perturbing them simultaneously. We are particularly interested in modulating the expression levels of genes that can be manipulated in DCs and thereby improving the immunogenic potency of DC vaccines (Figure 10A). For instance, one can upregulate the expression level of NF-κB through electroporating DCs with mRNA encoding constitutively active IKKβ. In addition, one could downregulate IκBα using microRNAs that repress gene expression at the post-transcriptional level to increase the level of NF-κB in DCs. One could also increase the levels of DC-secreted cytokines (such as IL-6 and IL-8) that are involved in the T-cell response through electroporation of the corresponding mRNAs into DCs.

As shown in Figure 10B left column, directly increasing the expression of free NF-κB shows dominant effects on the upregulation of memory T cells for any combination with other parameters. However, this manipulation is experimentally difficult, as NF-κB is a protein complex and requires the presence of subunits and heterodimerization to be functional. Alternatively, it is experimentally achievable by manipulating the expression of NF-κB regulators. Simultaneously modulating the expression of IKKβ (upregulation) and IκBα (downregulation) can result in an effective increase of memory T cells (Figure 10B middle column). Modulating the expression of IκBα has stronger effects than modulating IKKβ. Compared to single modulation of the NF-κB regulators, the combined modulation increases the population of memory T cells by 80–120%. Modulating the levels of IκBα together with DC-secreted cytokines (IL-6 and IL-8) is also an effective regulation of memory T cells but shows different dynamics compared to the modulation of both NF-κB regulators (i.e., IKKβ and IκBα) (Figure 10B middle column). The cytokines are less influential than IκBα in regulating the T-cell response, as modulating IκBα directly affects the expression of free NF-κB that can regulate the expression of multiple cytokines and membrane proteins (i.e., IL-6, IL-8, IL-12, and CD70) involved in T-cell activation. When we manipulated the expression level of IKKβ with a cytokine, the effects on increasing memory T cells are mild (Figure 10B right column). When both cytokines were simultaneously modulated, the effect on increasing memory T cells is also limited. This is due to the reason that T-cell activation also depends on other proteins (i.e., IL-12 and CD70), and their unchanged levels act as a limiting factor in the increase of memory T cells. This is in line with findings in a biological model system examining CTL-priming and memory formation in the presence or absence of T-cell help, where it was shown that direct cell-cell contact was crucial and soluble factors were not sufficient (Hoyer et al., 2014). Taken together, the simulations predicted that combined manipulation of IκBα and cytokines is an efficient strategy for increasing memory T cells in experiments, and such manipulation shows better performance than manipulating the expression levels of only the NF-κB regulators or the cytokines.