Interleukin-15 Signaling in HIF-1α Regulation in Natural Killer Cells, Insights Through Mathematical Models

Natural killer (NK) cells belong to the first line of host defense against infection and cancer. Cytokines, including interleukin-15 (IL-15), critically regulate NK cell activity, resulting in recognition and direct killing of transformed and infected target cells. NK cells have to adapt and respond in inflamed and often hypoxic areas. Cellular stabilization and accumulation of the transcription factor hypoxia-inducible factor-1α (HIF-1α) is a key mechanism of the cellular hypoxia response. At the same time, HIF-1α plays a critical role in both innate and adaptive immunity. While the HIF-1α hydroxylation and degradation pathway has been recently described with the help of mathematical methods, less is known concerning the mechanistic mathematical description of processes regulating the levels of HIF-1α mRNA and protein. In this work we combine mathematical modeling with experimental laboratory analysis and examine the dynamic relationship between HIF-1α mRNA, HIF-1α protein, and IL-15-mediated upstream signaling events in NK cells from human blood. We propose a system of non-linear ordinary differential equations with positive and negative feedback loops for describing the complex interplay of HIF-1α regulators. The experimental design is optimized with the help of mathematical methods, and numerical optimization techniques yield reliable parameter estimates. The mathematical model allows for the investigation and prediction of HIF-1α stabilization under different inflammatory conditions and provides a better understanding of mechanisms mediating cellular enrichment of HIF-1α. Thanks to the combination of in vitro experimental data and in silico predictions we identified the mammalian target of rapamycin (mTOR), the nuclear factor-κB (NF-κB), and the signal transducer and activator of transcription 3 (STAT3) as central regulators of HIF-1α accumulation. We hypothesize that the regulatory pathway proposed here for NK cells can be extended to other types of immune cells. Understanding the molecular mechanisms involved in the dynamic regulation of the HIF-1α pathway in immune cells is of central importance to the immune cell function and could be a promising strategy in the design of treatments for human inflammatory diseases and cancer.

1. Introduction

As effector lymphocytes of innate immunity, natural killer (NK) cells are involved in the host defense against microbial infections and cancer (1). Sensing their environment, NK cells respond to cellular alterations including those caused by infections, cellular stress, and transformation (2).

Interleukin-15 (IL-15), produced by monocytes, macrophages and dendritic cells, critically regulates NK cell survival and activation (3, 4). While expression of IL-15 is low under homeostatic conditions, it is upregulated in inflammation (5). Upon receptor binding, IL-15 initiates Janus kinase/signal transducer and activator of transcription (JAK/STAT) signaling. This promotes growth of NK cells and enhances their ability to respond to activation. Activated NK cells infiltrate tissues containing pathogen-infected or malignant cells, resulting in their recognition and direct killing (6–8).

Sites of infection or cellular transformations are often characterized by inflammatory hypoxia. Thus, NK cells must adapt and respond under conditions of low oxygen tension. The critical cellular dependence of survival on oxygen led to the early evolution of adaptive cellular responses to hypoxia. Cellular adaptation to hypoxia is primarily orchestrated by the hypoxia inducible factor (HIF) family of transcription factors (9). To date, three HIF family members have been identified (HIF-1, HIF-2, and HIF-3) of which HIF-1 is the best characterized (10). Two subunits, HIF-1α and HIF-1β, form the transcriptionally active HIF-1 complex. The α-subunit is post-translationally hydroxylated by oxygen-sensitive prolyl hydroxylases (PHDs) which mark the protein for ubiquitination and continuous proteasomal degradation. A decrease in cellular oxygen availability stabilizes HIF-1α allowing its dimerization with HIF-1β. The dimer translocates to the nucleus, binds to hypoxia-response elements in promoters of adaptive genes, and activates their expression.

In immune cells, including T lymphocytes and myeloid cells, cellular activation of HIF-1α has also been reported to occur in an oxygen-independent manner during inflammation triggered by infection and cancer, and to involve transcriptional in addition to post-translational mechanisms (11). At sites of tissue damage and infection, both inflammation and decreased oxygen availability result in HIF-1α stabilization and its nuclear translocation. Recent insights from models of solid tumors in mice with an NK cell specific knockout of the HIF-1α gene and from chemical inhibition of HIF-1α in human NK cells (12, 13) suggest that HIF-1α limits NK cell anti-tumor activity.

In the past 15 years, several mathematical models for HIF-1α regulation based on systems of ordinary differential equations (ODEs) have been proposed (14–19). A review up to 2013 is given in (20). Nguyen and coauthors (19) have investigated the dynamics of the HIF-1α pathway, combining a mathematical mechanistic model and experimental analysis for human embryonic kidney 293 (HEK-293) cells. Their model studies accumulation of HIF-1α in hypoxia and its degradation in normoxia, considering hydroxylation of HIF-1α mediated through both prolyl hydroxylases and asparaginyl hydroxylase FIH (factor inhibiting HIF). Fábián et al. (10) have highlighted the importance of using system biology and mathematical modeling for understanding HIF signaling. Although lacking the comparison with experimental data, Fábián's models allowed to test different hypotheses on the HIF network, concluding that the negative feedback induced by PHDs plays a major role in triggering oscillations in the HIF-1α dynamics.

This work combines mathematical modeling and experimental analysis to understand processes regulating the levels of HIF-1α mRNA and protein in NK cells. The proposed mathematical model considers key features of HIF-1α regulation and is formulated as a system of non-linear ODEs with positive and negative feedbacks. In our in vitro studies, we isolated human peripheral NK cells and studied their behavior simulating hypoxic and inflammatory conditions, which were produced by the hypoxia-mimicking agent dimethyl-oxalyl glycine (DMOG) and the pro-inflammatory cytokine IL-15, respectively. Experimental trials were designed to collect time series data of HIF-1α protein expression and its upstream regulators in order to calibrate the mathematical model. Parameter estimation was performed by means of numerical methods based on a multiple shooting approach for dynamic systems and a generalized Gauss-Newton method for optimization. Our approach does not only explain experimental observations on HIF-1α dynamics but also allows to ask questions and test hypotheses with the help of in silico experiments. For example, we investigated how HIF-1α levels depend on the regulation of other upstream proteins, and identified the signal transducer and activator of transcription 3 (STAT3), the mammalian target of rapamycin (mTOR) and the nuclear factor-κB (NF-κB) as critical regulators. Further, we studied HIF-1α stabilization in dependence of DMOG-mediated PHD/FIH inhibition, determining a non-linear relation between HIF-1α levels and DMOG concentration. Our model provides new insights into the mechanisms mediating accumulation of HIF-1α in NK cells, by (i) highlighting the synergistic effects of IL-15 and chemical hypoxia, and (ii) suggesting that NF-κB and STAT3 are fundamental regulators of IL-15 induced HIF-1α enrichment.

2. Materials and Methods

2.1. NK Cell Purification and Cell Culture

The study was reviewed and approved by the Medical Ethics Commission II of the Medical Faculty Mannheim, Heidelberg University (2014-500N-MA). NK cells were isolated from buffy coats obtained through the local Red Cross Blood Donor Service (NK-Cell Isolation Kit, Miltenyi Biotec GmbH, Bergisch Gladbach, Germany). The purity of NK cells was determined by flow cytometry.

Freshly isolated NK cell preparations with a phenotype of ≥95% CD56⁺CD3⁻ and ≤1% each CD3⁺, CD14⁺, CD15⁺, and CD19⁺ were judged as pure and were further cultivated as previously described (21). In brief, cells were plated at a density of 10⁶ cells/mL in RPMI 1640 medium (Sigma-Aldrich Chemie GmbH, Merck KGaA, Darmstadt, Germany) supplemented with 10% fetal bovine serum (FBS) and 2 mM L-glutamine and maintained in a standard tissue culture incubator (37°C, 5% CO₂, 21% O₂, normoxia, standard condition). The cell permeable pan-hydroxylase inhibitor DMOG (Selleck Chemicals, Houston, TX, USA) was used to mimic hypoxia. The viability of the cells was determined by tryptan blue staining and was ≥95% (Countess, Invitrogen, ThermoFisher, Waltham, MA, USA).

2.2. In vitro Treatments

Freshly isolated NK cells were maintained overnight under standard conditions and were stimulated with human recombinant IL-15 (45 ng/mL, PeproTech, NJ, USA), DMOG (20 μM, Selleck Chemicals), rapamycin (25 nM, Merck Chemicals GmbH, Darmstadt, Germany), STAT3 inhibitor (S3I-201, 200 μM, Merck Chemicals GmbH), or DMSO (Sigma-Aldrich Chemie GmbH) as control, on the next day for the indicated time periods. Protein concentrations in cell lysates were determined on a Direct Detect^® infrared spectrometer (Merck Millipore) according to the manufacturer's instructions.

2.3. Western Blotting

Total cell extracts were prepared by resuspending 3 × 10⁶ NK cells in 100 μL NP-40 lysis buffer (50 mM Tris-HCl, pH 7.5, 120 mM NaCl, 20 mM NaF, 1 mM EDTA, 6 mM EGTA, 15 mM sodium pyrophosphate, 1 mM PMSF, 0.1% Nonident P-40). Fifteen minutes of cell lysis on ice was followed by centrifugation for 20 min at 14,000 × g. Cleared lysates were analyzed directly by SDS-PAGE and Western Blotting. Briefly, equal amounts of protein were separated by SDS-PAGE, transferred to nitrocellulose membranes (Thermo Fisher), blocked in 5% dry milk powder dissolved in 1×PBS-T, and then probed with primary antibody and HRP-conjugated secondary antibody (Santa Cruz Biotechnology, Dallas, TX, USA). Proteins were visualized using Enhanced Chemiluminescent solution (Thermo Fisher) and FUSION Vilber imager (Eberhardzell, Germany). The intensity of signals was quantified by densitometric analysis using the image analysis software ImageJ (Version 1.51j8). The value for HIF-1α was normalized to that for β-Actin. Anti-HIF-1α (# 2185) was obtained from Abcam (Cambridge, UK) and Anti-Actin (8H10D10) from Cell Signaling Technology (Frankfurt am Main, Germany). Representative experiments out of three performed are shown.

2.4. MILLIPLEX Immunoassay

The MILLIPLEX MAP Multi-Pathway Signaling Phosphoprotein Kit 48-680MAG was used according to the manufacturer's protocol (Merck Millipore). Total NK cell extracts were diluted with MILLIPLEX MAP Assay Buffer to reach the protein concentration of 10 μg of total protein/well. Mixed magnetic beads were added to each well. To appropriate wells, 25 μL of Assay Buffer (background control), 25 μL of NK cell sample lysates and 25 μL of control cell lysates were added in duplicates. The plate was sealed and incubated overnight (20 h) at 4°C on a plate shaker (750 rpm). After incubation and washing, 25 μL of Detection Antibody were added to each well. This incubation step was followed by addition of 25 μL Streptavidin-Phycoerythrin and incubation on the shaker. After resuspending the beads in 150 μL of Assay Buffer, the plate was read on a MAGPIX system (Luminex). Signals of phosphorylation of STAT3 and AKT were expressed as background-corrected median fluorescence intensities.

2.5. Modeling the Regulatory Network

The mathematical approach used in this study is based on a system of autonomous non-linear ODE, which can be in general written as

$\begin{array}{ll}\{\begin{array}{ccc}{y}^{\prime }(t)& =f(y,u,p),& t\in \left[t_0,t_f\right]\subset ℝ\\ y(t_0)& =y_0(p).& \end{array}& (1)\end{array}$

with states y, controls u and parameters p. The vector $y(t)\in {ℝ}^{n_y}$ indicates the “state of the system,” that is, the concentration of the considered proteins, complexes and mRNA at time t ∈ [t₀, t_f] ⊂ ℝ, and $y_0\in {ℝ}^{n_y}$ is the initial state. The parameter vector $p\in {ℝ}^{n_p}$ contains non-negative constants describing the biochemical reaction rates (such as production, degradation, binding, etc.) in the system. The time-dependent experimental controls $u(t)\in {ℝ}^{n_u}$ represent cell treatments, specifically with IL-15, DMOG, or other protein inhibitors. In this work we do not discuss basic theoretical properties of the solutions to system 1, such as existence and uniqueness of a global solution, or invariance of the positive cone of ${ℝ}^{n_y}$. All these properties can be proven by applying elementary results and methods in ODE theory [see e.g., (22)], and we assume them to hold true in this manuscript. In the following we explain in detail the model assumptions and the resulting equations. Model variables and parameters are given in Table S1, and Tables 1, 2, respectively. A diagram of the regulatory network is shown in Figure 1.

TABLE 1

Table 1. Parameter description and values used for the mathematical model (2)–(11).

TABLE 2

Table 2. Parameter description and values used for the mathematical model (2)–(11).

FIGURE 1

Figure 1. Diagram for HIF-1α regulatory network in NK cells, corresponding to model Equations (2)–(11). We study the interplay of HIF-1α, IL-15, mTOR, NF-κB, and STAT3 in normoxia and hypoxia. A signaling cascade starting with IL-15 activates NF-κB, STAT3, and AKT. This in turn activates mTOR, influencing HIF-1α mRNA and HIF-1α protein levels. In normoxia, HIF-1α is hydroxylated via FIH and PHD. Here, we consider only the FIH-mediated asparaginyl-hydroxylated HIF-1α (HIF-1α-aOH) and assume that a small fraction of HIF-1α-aOH can be dehydroxylated. In normoxia, hydroxylated HIF-1α is degraded via Von Hippel-Landau protein (not considered in the model). In hypoxia, HIF-1α accumulates and binds to HIF-1β, building the HIF-1 complex. The latter is responsible for both a positive and a negative feedback on HIF-1α, via activation of NF-κB and via inhibition of mTOR and upregulation of PHD. Black arrows indicate protein activation, translation or formation/dissociation of protein complexes. For simplicity, basal degradation of molecules is not depicted here. The blind-ended arrow indicates inhibition. Blue arrows indicate external regulation due to further stimuli not specifically considered in the mathematical model. Yellow arrows indicate HIF-1α hydroxylation via FIH or PHD.

Recent results (24, 25) showed the connection between IL-15 and mTOR activity in NK cells, indicating that the AKT-mTOR pathway is indispensable for efficient cell activity and immune functions of NK cells. We therefore focused on the latter signaling pathway, neglecting other cascades [such as Ras-Raf-MEK and JAK/STAT5 (26)] which are also known to be initiated by IL-15. Further, IL-15 stimulation in neutrophils and human peripheral blood lymphocytes has been shown to activate NF-κB and STAT3 (27–29). All in all, we assumed that IL-15 (y₁) activates AKT (y₂), NF-κB (y₇), and STAT3 (y₈). For a general formulation we further assumed that IL-15 enters the system at constant rate a₁ and decays at rate d₁, that is, the IL-15 dynamics is given by

$\begin{array}{ll}{y}_1^{\prime }(t)=a_1-d_1{y}_1.& (2)\end{array}$

Activation of AKT is assumed to occur via IL-15 (30) (activation rate k₁), other external mediators (basal activation rate a₂), and also via STAT3 (maximal rate k_S) (31). We further assumed that AKT constantly decays at rate d₂, yielding

$\begin{array}{ll}{y}_2^{\prime }(t)=a_2+k_1{y}_1+k_S\frac{y_8^{n_2}}{{\xi }_{28}^{n_2}+y_8^{n_2}}-d_2{y}_2.& (3)\end{array}$

Following the literature (32) we assumed that AKT activates mTOR (y₃) at rate k₂. mTOR basal activation and decay rate are denoted by a₃ and d₃, respectively. The inhibitory effect that hypoxia has on mTOR (33) is included in the model by means of a negative feedback regulated by the HIF-1 complex (y₆). Hence, for the mTOR dynamics we obtained

$\begin{array}{ll}{y}_3^{\prime }(t)=\left(a_3+k_2{y}_2\right)\frac{{\alpha }_1}{{\alpha }_2+y_6}-d_3{y}_3.& (4)\end{array}$

We denote phosphorylated STAT3 by y₈. STAT3 basal activation and decay rates are a₈ and d₈, respectively. Both IL-15 and mTOR are known to induce phosphorylation of STAT3 (29, 34), here assumed to occur at rate k₆ and k₈, respectively. All in all, the differential equation for STAT3 reads

$\begin{array}{ll}{y}_8^{\prime }(t)=a_8+k_8{y}_3+k_6{y}_1-d_8{y}_8.& (5)\end{array}$

The last protein upstream of HIF-1α that we considered in our model is NF-κB, for which we assumed basal activation (at rate a₇) and decay (d₇). NF-κB is further activated via IL-15 (27, 28) (activation rate k₇), via mTOR (35) (k₁₅) and via the HIF-1 complex (36, 37) (k₁₄), yielding

$\begin{array}{ll}{y}_7^{\prime }(t)=a_7+k_7{y}_1+k_{14}{y}_6+k_{15}{y}_3-d_7{y}_7.& (6)\end{array}$

HIF-1α mRNA basal synthesis and degradation are defined to occur at rate a₉ and d₉, respectively. Further, we assumed that HIF-1α mRNA is regulated by NF-κB (at rate k₉) and STAT3 (at rate k₃),

$\begin{array}{ll}{y}_9^{\prime }(t)=a_9+k_9{y}_7+k_3{y}_8-d_9{y}_9.& (7)\end{array}$

Following previous results (19), we assumed that asparaginyl hydroxylase FIH is at steady state (φ), whereas PHDs are upregulated by HIF-1 complex and we approximated their dynamics with quasi-steady state assumptions (see section S1 for detailed explanation). Further, we assumed that HIF-1α mRNA is translated at rate k_α and HIF-1α protein decays at rate d₄. We denote by K_O2 the oxygen-dependent binding force of FIH/PHD and HIF-1α (cf. section S1). In normoxia, HIF-1α is hydroxylated via FIH (assumed at maximal rate k₁₀) and via PHD (maximal rate k₁₃). The dynamics of HIF-1α protein (y₄) is thus given by

$\begin{array}{lll}{y}_4^{\prime }(t)& =& k_{\alpha }{y}_9-d_4{y}_4-k_4{y}_4{y}_5+k_5{y}_6-k_{13}{K_{\text{O}}}_2(\Delta y_6+a_{11})\frac{y_4}{{\xi }_{44}+y_4}\\ \hfill & & -k_{10}{K_{\text{O}}}_2\phi \frac{y_4}{{\xi }_4+y_4}+k_{11}{y}_{10}.\hfill & (8)\end{array}$

In accordance with previous studies (19, 38), we assumed that asparaginyl-hydroxylated HIF-1α (HIF-1α-aOH, here denoted by y₁₀) can be subsequently hydroxylated via PHD and then degraded, whereas prolyl-hydroxylated HIF-1α (HIF-1α-pOH) is quickly degraded. We henceforth neglected the dynamics of HIF-1α-pOH. Further, we assumed that there is some probability for HIF-1α-aOH dehydroxylation [cf. (19)]. The resulting dynamics of HIF-1α-aOH is given by

$\begin{array}{lll}{y}_{10}^{\prime }(t)& =& k_{10}{K_{\text{O}}}_2\phi \frac{y_4}{{\xi }_4+y_4}-k_{12}{K_{\text{O}}}_2(\Delta y_6+a_{11})\frac{y_{10}}{{\xi }_{10}+y_{10}}\\ & & -k_{11}{y}_{10}-d_{10}{y}_{10}.& (9)\end{array}$

HIF-1β (y₅) is constitutively expressed by the cells (synthesis rate a₅), independently of the oxygen conditions (10). In hypoxia, HIF-1α accumulates and binds to HIF-1β (at rate k₄) forming the transcriptional complex HIF-1, which can dissociate (rate k₅). Hence, for HIF-1β and the HIF-1 complex we obtained the equations,

$\begin{array}{ll}{y}_5^{\prime }(t)=a_5-k_4{y}_4{y}_5+k_5{y}_6-d_5{y}_5,& (10)\end{array}$

and

$\begin{array}{ll}{y}_6^{\prime }(t)=k_4{y}_4{y}_5-k_5{y}_6-d_6{y}_6,& (11)\end{array}$

respectively. For model calibration and comparison with collected experimental data we extended the model (2)–(11) to include DMOG or other protein inhibitors. Details are given in the Supplementary Material (section S2).

2.6. Numerical Simulations

For the numerical integration of the non-linear ODE system (2)–(11) and numerical simulations shown in this manuscript we used the Runge-Kutta formula (4,5) and (3,2) in MATLAB^® version 9.4 [routines ode45 and ode23s (39)], as well as a multistep Backward Differentiation Formula (BDF) method with variable step size and order control. The latter was implemented by mean of the solver DAESOL (40, 41) in VPLAN (42), a software for simulation, parameter estimation and optimum experimental design for non-linear processes described by differential equations.

We assumed that the initial load of IL-15 in primed cells is y₁(0) = 1, whereas for unstimulated cells y₁(0) = 0. Collected experimental data for HIF-1α, STAT3 and AKT were normalized with respect to measurements in untreated cells at the beginning of each experiment and used for model calibration. In setting the initial conditions, we normalized AKT, mTOR, NF-κB, STAT3, HIF-1α mRNA and HIF-1β with respect to the concentration at the beginning of each in silico experiment, meaning that y_j(0) = 1, for j = 2, 3, 5, 7, 8, 9. The total HIF-1α level was normalized with respect to the initial measurement, corresponding to untreated cells in normoxia, hence we set y₄(0) + y₆(0) + y₁₀(0) = 1. As cells were pre-cultivated in normoxia, we assumed that at the beginning of our observations (t = 0 h) most of HIF-1α is hydroxylated, hence, we set y₄(0) = 0.05, y₆(0) = 0.05 and y₁₀(0) = 0.9.

2.7. Model Calibration

2.7.1. Comparison With Experimental Data

For comparison with experimental data, the solution y(·) = y(·|u, p) of the mathematical model (1) is associated with an m-vector of observables,

$\begin{array}{l}g(t,y(t),u,p).\end{array}$

Typically it is not possible to observe all states, and in general g(·) is a non-linear function de-pending on the states and parameters. Given an experimental setting u^ex(·), for each observable $g_i(t_j^{ex},y^{ex}(t_j^{ex}),u^{ex},p),i=1,\dots ,m$, at measurement times $t_j^{ex}\in \left[t_0,t_f\right],j=1,\dots ,k^{ex}$, experimental data ${\eta }_{ij}^{ex}$ are collected in each experiment ex. Experimental measurements contain additive noise

$\begin{array}{ll}{\eta }_{ij}^{ex}=g_i(t_j^{ex},y^{ex}(t_j^{ex}),u^{ex},p)+{\epsilon }_{ij}^{ex},& (12)\end{array}$

where the errors ${\epsilon }_{ij}^{ex}$ are assumed to be statistically independent and normally distributed with zero mean value and variances ${({\sigma }_{ij}^{ex})}^2$.

The experimental settings considered in this study present “perturbation” experiments, ex = 1, …, n_e, which allow to investigate perturbations of the cellular processes from their equilibrium conditions. Each perturbation experiment corresponds to a different control u = u^ex in the mathematical model (1), and the same quantities are measured in each experiment. It is biologically reasonable to assume that an unperturbed system is at the steady state. This corresponds in our case to unstimulated NK cells (u^ex = u⁰), hence to the initial condition y₀. The steady state can be mathematically determined considering the dynamics of untreated cells and setting this in equilibrium (43, 44). Hence, the initial condition of the system y₀ and the model parameters must satisfy the steady state constraint,

$\begin{array}{ll}0=f(y_0,u^0,p).& (13)\end{array}$

2.7.2. Parameter Estimation Problem

In general, given a mathematical model (1) and experimental measurements, the goal of model calibration is to determine the model parameters in (1) from the collected data. This reduces to an optimization problem of minimizing the discrepancy between model observables and the experimental data using a particular metric (45). The weighted least squares functional is known to deliver a maximum likelihood estimate for the unknown parameters (46, 47). For the calibration of the parameters of the regulatory network (2)–(11) we used a weighted l₂-norm of the measurement errors,

$\begin{array}{ll}{l}_2(y,p):=\frac{1}{2}{\displaystyle \sum _{i,j,ex}}\frac{{({\eta }_{ij}^{ex}-g_i(t_j^{ex},y^{ex}(t_j^{ex}),u^{ex},p))}^2}{{({\sigma }_{ij}^{ex})}^2},& (14)\end{array}$

and further included a priori information p₀ by adding a Tikhonov regularization term (48)

$\begin{array}{l}L(p,p^0,\lambda ):=\frac{1}{2}{\displaystyle \sum _{\tilde{m}=1}^{n_p}}\frac{{(p_{\tilde{m}}-p_{\tilde{m}}^0)}^2}{{\lambda }_{\tilde{m}}^2},\end{array}$

where the vector $\lambda \in {ℝ}^{n_p}$ controls the amount of regularization per parameter. Moreover we incorporated additional information about parameters and states (initial conditions, steady states, etc.) in the parameter estimation problem by formulating equality constraints (49, 50). We estimated parameters by solving the following multiple experiment parameter estimation (PEP) problem

$\begin{array}{l}\text{(PEP) }\{\begin{array}{c}{\displaystyle \underset{y(·),p}{\text{min}}}\frac{1}{2}{\displaystyle {\sum }_{i,j,ex}}\frac{{({\eta }_{ij}^{ex}-g_i(t_j^{ex},y^{ex}(t_j^{ex}),u^{ex},p))}^2}{{({\sigma }_{ij}^{ex})}^2}+\frac{1}{2}{\displaystyle {\sum }_{\tilde{m}=1}^{n_p}}\frac{{(p_{\tilde{m}}-p_{\tilde{m}}^0)}^2}{{\lambda }_{\tilde{m}}^2},\hfill \\ \text{s.t. }{\left(y^{ex}\right)}^{\prime }(t)=f(y^{ex},u^{ex},p),t\in \left[t_0,t_f\right],y^{ex}(t_0)=y_0^{ex}(p),\hfill \\ r^{ex}(y^{ex}(t_1^{ex}),y^{ex}(t_2^{ex}),\dots ,y^{ex}(t_k^{ex}),p)=0,ex=1,\dots ,n_e.\hfill \end{array}\end{array}$

2.7.3. Numerical Methods for Parameter Estimation

For least squares minimization, as those in (PEP), a frequently adopted approach is the derivative-based iterative Gauss-Newton method (45, 51). In this work we applied an “all-at-once” parameter estimation method based on a direct multiple shooting approach for dynamic systems (51) and a generalized Gauss-Newton method for optimization (50, 51). This is a boundary value problem approach, in which the system of differential Equations (1) is discretized including boundary conditions. The discretized system is treated as a non-linear constraint of the least squares objective function (52). For the multiple shooting approach, a suitable grid of multiple shooting nodes was chosen and, at each multiple shooting grid point, the values of the state variables were included as additional optimization variables. On each subinterval an additional initial value problem was solved. To maintain the continuity and feasibility of the solution, we included additional matching conditions (50, 52). The splitting of the integration interval leads to a numerically stable system (52). The resulting finite dimensional non-linear constrained least squares problem can be formally written as

$\begin{array}{ll}{\displaystyle \underset{s,p}{\text{min}}}\frac{1}{2}\Vert F_1(s,p){\Vert }_2^2,\text{ s.t. }{F}_2(s,p)=0,& (15)\end{array}$

where the constraints F₂(s, p) = 0 include the multiple shooting parameterization of the dynamical model and s denotes the vector of states in the parametrized model. The problem (15) was solved by a generalized Gauss-Newton method. At each iteration of the Gauss Newton method

$s^{k+1}=s^k+t^k\Delta s^k,p^{k+1}=p^k+t^k\Delta p^k,t^k\in ]0,1]$

the increments Δs^k, Δp^k solve the linearized problem

$\begin{array}{ll}\{\begin{array}{c}{\displaystyle \underset{\Delta s,\Delta p}{min}}\frac{1}{2}\Vert F_1(s^k,p^k)+\frac{\partial F_1(s^k,p^k)}{\partial s}\Delta s+\frac{\partial F_1(s^k,p^k)}{\partial p}\Delta p{\Vert }_2^2,\\ \text{ s.t. }{F}_2(s^k,p^k)+\frac{\partial F_2(s^k,p^k)}{\partial s}\Delta s+\frac{\partial F_2(s^k,p^k)}{\partial p}\Delta p=0.\end{array}& (16)\end{array}$

For the case considered in this work, the linearized problem (16) shows special structures due to multiple experiments and multiple shooting approaches. These structures are efficiently exploited in a tailored linear algebra method for the solution of (16). A numerical analysis of the well-posedness of the problem and an assessment of the error of the resulting parameter estimates were performed at the solution of the problem (PEP), based on the analysis of the corresponding sensitivity (or Jacobian) matrix J

$\begin{array}{ll}J=J(s,p)=\left(\begin{array}{c}{J}_1(s,p)\\ J_2(s,p)\end{array}\right)=\left(\begin{array}{cc}\frac{\partial F_1(s,p)}{\partial s}& \frac{\partial F_1(s,p)}{\partial p}\\ \frac{\partial F_2(s,p)}{\partial s}& \frac{\partial F_2(s,p)}{\partial p}\end{array}\right).& (17)\end{array}$

Further details can be found in (49, 50). In particular, we computed the linear approximation of the variance-covariance matrix for the constrained parameter estimation problem (PEP)

$\begin{array}{ll}Cov=Cov(s,p)=J^{+T}\left(\begin{array}{cc}\text{diag}{\left\{{\sigma }_{ij}^{ex}\right\}}_{i,j,ex}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right)J^{+},& (18)\end{array}$

and the standard deviations of states and parameters as square root of the corresponding diagonal elements of the matrix Cov. In (18) the matrix J⁺ denotes the generalized inverse of the sensitivity matrix J(s, p), that is J⁺JJ⁺ = J, and 0 denotes the zero matrices of the corresponding dimensions.

The above sketched methods are implemented in the software VPLAN (42), which includes the parameter estimation software PARFIT (53, 54). VPLAN was used for parameter estimation in this study.

2.7.4. Initial Guesses

Before starting the actual optimization procedure (PEP), we determined initial guesses in a reasonable scale for the parameters that needed to be estimated. Good initial guesses are important for convergence of the parameter estimation methods used in this work. In certain cases prior information p⁰ for initial guesses is found in the literature (cf. Table 2). When results published in previous studies did not help, we applied homotopy related methods (55) and pre-estimated the parameters using Wolfram Mathematica^® version 11.3 with the Ndsolve and manipulate routines, as well as DAESOL in VPLAN.

2.7.5. Sensitivity Analysis

In a second step, before starting parameter estimation, we performed a local sensitivity analysis at parameter values p⁰ and corresponding solutions y^ex(·|u^ex, p⁰), ex = 1, …, n_e of ODE systems. The goal of the local sensitivity analysis is to find those parameters, which can be estimated reliably with the model and the given measurements. The local sensitivity analysis is performed using the sensitivity matrix J⁰ = J(s⁰, p⁰). If the sensitivity matrix J⁰ has almost linear dependent columns, then it is ill-conditioned and the parameter vector is badly identifiable (48, 56). We computed the singular value decomposition of the sensitivity matrix J⁰. The reciprocal of the minimal singular value yields the collinearity index, which, if it is large, indicates that the sensitivity matrix has almost linear dependent columns (48, 56, 57). In this work we set a minimal threshold 0.1 for rejecting small singular values and for determining the subset of parameters corresponding to almost linear independent columns of J⁰. This allowed to fix 13 parameters which correspond to almost linear dependent columns of J⁰. The remaining 25 parameters are reliably identifiable and were estimated by mean of (PEP).

2.8. Model Robustness

We implemented extrinsic and intrinsic stochastic perturbations using a Monte Carlo analysis. For extrinsic perturbations we varied the input stimulation IL-15 (y₁(t)) ±25% around its original value. We measured the effect of these perturbations on total HIF-1α in the model accounting for DMOG treatment and IL-15 stimulation at the steady-state level. For intrinsic perturbations we varied specific parameters (a₂, a₃) and measured the effect of these perturbations on total HIF-1α in the model accounting for DMOG treatment and IL-15 stimulation at the steady state level. We implemented the intrinsic and extrinsic stochastic perturbations by varying the specific elements 25% around their original value and sampling 1000 times from a uniform distribution. The Monte Carlo analysis was implemented in Wolfram Mathematica^® version 11.2.

3. Results

3.1. IL-15-induced HIF-1α Protein Accumulation in Peripheral NK Cells

In cells exposed to hypoxia, the stabilization and activation of HIF-1α is well characterized (58). Instead of manipulating the oxygen tension to induce HIF-1α, pharmacological inhibitors of HIF-1α hydroxylation can be used as well. We used a cell-permeable pan-hydroxylase inhibitor, DMOG, to inhibit oxygen-sensitive hydroxylases that target HIF-1α for proteasomal degradation and its transcriptional inactivation. After preincubation under normoxia for 16 h, peripheral NK cells were stimulated with the pro-inflammatory cytokine IL-15 in the presence of DMOG for different time periods. Whole cell extracts were prepared and the response of NK cells to IL-15 stimulation and DMOG treatment was assessed by evaluating the expression of HIF-1α measured by Western Blot analysis (Figure 2). The expression of β-actin was monitored to confirm equal loading.

FIGURE 2

Figure 2. In the presence of DMOG, stimulation of NK cells with IL-15 induces the accumulation of HIF-1α protein. Whole cell extracts were used for immunoblotting of HIF-1α. Equal gel loading was confirmed by β-actin expression. (A) After preincubation under normoxia for 16 h, NK cells were treated with IL-15 and DMOG for indicated additional time periods. (B) In addition to IL-15 and DMOG, NK cells were treated with rapamycin (Rapa), which reduces HIF-1α protein accumulation.

HIF-1α expression was barely detectable in the first 3 h of IL-15 and DMOG stimulation. However, after 4 h we detected accumulation of HIF-1α protein, which further increased over 8 h and was maintained to at least 27 h, although to a lesser extent than at 8 h (Figure 2A).

IL-15 signaling in NK cells through the kinase mTOR has previously been reported to be essential for their expansion in the bone marrow and sustained activation (24). Moreover, among other signaling pathways, the PI3K/mTOR pathway has been linked to the induction of HIF-1α protein expression in immune cells, including T lymphocytes (59). To study the role of mTOR in HIF-1α protein expression in NK cells, we stimulated the cells with IL-15 and DMOG in the presence or absence of pharmacological mTOR inhibitor rapamycin. As shown in Figure 2B, mTOR inhibition reduced HIF-1α levels. Nevertheless, HIF-1α signals remained detectable, pointing to other upstream regulators of HIF-1α protein accumulation in IL-15 stimulated NK cells.

3.2. Model Parameters

The model (2)–(11) has been calibrated on time series for AKT, STAT3 and HIF-1α collected from NK cells under different experimental conditions (Table S2), namely under stimulation/treatment with: (i) IL-15; (ii) DMOG; (iii) IL-15 + DMOG + rapamycin; (iv) DMOG + IL-15 + S3I-201. We interpret (i)–(iv) as “perturbation experiments” from the initial (equilibrium) condition. Biological sense and previous literature on mathematical modeling of cellular dynamics (43, 44) suggest that it is important to assume that untreated cells are in the steady state. Such an assumption yields 10 algebraic equations (steady state constrains) on the parameter, which together with the collected time series data of the four perturbation experiments (39 data points, cf. Table S2) were used to calibrate the model. Table 2 reports the 25 estimated parameters and Figure 3A shows the model fit. Model parameters either estimated or fixed from previous literature are reported in Tables 1, 2. With these parameter values we ran an in-silico experiment for NK cells stimulated with IL-15 and treated with DMOG. The obtained numerical simulations were used to validate our model, comparing the model predictions with collected experimental data for AKT, STAT3 and HIF-1α time series of NK cells stimulated with IL-15 and treated with DMOG (Figure 3B). To depict the statistical significance of the parameter estimates, Table 2 also reports the standard deviation for each estimated parameter value.

FIGURE 3

Figure 3. The mathematical model (2)–(11) explains collected time series for HIF-1α, STAT3, and AKT. Parameters used for numerical simulations are given in Tables 1, 2. (A) Model calibration results: comparison of numerical simulations (continuous curves) and collected experimental data (dots ± S.E.) of total HIF-1α (y₄ + y₆ + y₁₀, red curves/dots), STAT3 (blue) and AKT (black). The model is fitted to data collected in different experimental settings: (upper left) IL-15-stimulated NK cells; (upper right) DMOG treated cells; (lower left) IL-15-stimulated NK cells treated with DMOG and STAT3 inhibitor (S3I-201); (lower right) IL-15-stimulated cells treated with DMOG and with mTOR inhibitor rapamycin (Rapa). (B) Model validation results: comparison of numerical simulations (continuous curves) and collected experimental data (dots ± S.E.) of total HIF-1α (y₄ + y₆ + y₁₀, red curves/dots), STAT3 (blue), and AKT (black) for IL-15-stimulated cells treated with DMOG. Experimentally collected data points are reported in Table S2.

3.3. The Mathematical Model Explains the Dynamics of HIF-1α Accumulation

Besides the known hypoxia-induced HIF-1α stabilization and AKT-mTOR-mediated increase in protein translation, HIF-1α can also be induced through increased transcription involving activated transcription factors, among others STAT3 as shown in T lymphocytes (60) and B lymphocytes (61). To determine the role of AKT, mTOR and STAT3 as mediators of HIF-1α accumulation downstream of IL-15, we collected time series data (Figure 3, Table S2) for the phosphorylation status of AKT (Ser473) and STAT3 (Ser727), representing the activated forms of the proteins (62, 63). For optimal parameter estimation, we collected data from NK cells isolated from blood of the same donors, cultured in normoxia, chemical hypoxia (DMOG) and treated with a STAT3 or mTOR inhibitor.

The model is able to reproduce data collected for HIF-1α, STAT3 and AKT in different experimental settings (Figure 3). In particular, model predictions match time series data of HIF-1α protein expression and indicate that simultaneous exposure of NK cells to IL-15 and DMOG (Figure 3B) increases the levels of total HIF-1α, compared to HIF-1α levels in cells either stimulated with IL-15 or treated with DMOG (Figure 3A). Moreover, inhibition of mTOR or STAT3 leads to reduction of HIF-1α levels, suggesting that both proteins are involved in the regulation of IL-15 induced HIF-1α accumulation in DMOG treated cells. Model assumptions and calibration results (cf. Tables 1, 2) indicate that the external regulation of IL-15, NF-κB, STAT3, and HIF-1α-mRNA is negligible in order to explain collected time series in NK cells under the proposed experimental setting. Further, parameter estimation and model discrimination (results not shown here) suggest that DMOG reduces IL-15-mediated STAT3 activation (see section S2).

The collected data (cf. Table S2) further show that IL-15, DMOG or inhibition of either mTOR or STAT3 does not affect the levels of phosphorylated AKT (Ser473) in human NK cells. After preliminary steps in the parameter estimation procedure, we obtained $k_1\approx 10^{-5},k_S\approx 10^{-4}$ (cf. Table 1). The sensitivity analysis which followed indicated that the two values have small effects on the objective function in (PEP). This suggests that, the order of magnitude of k₁ and k_S being much smaller than those of all other parameters, the two parameters could be set to zero without much affecting the simulation results, and the AKT dynamics in Equation (3) can be simplified and described by a linear equation,

$\begin{array}{l}{y}_2^{\prime }(t)=a_2-d_2{y}_2.\end{array}$

With parameter values as indicated in Tables 1, 2 we ran numerical simulations of model (2)–(11) for NK cells under different experimental conditions. Figures 4, 5 show the simulation results for all¹ model components in normoxia without (N, O₂ = 21%) and with DMOG (20 μM, chemical hypoxia), treated with one (A panels) or two inhibitors at the same time (B panels). The level of total HIF-1α, which we define as the sum of HIF-1α protein, HIF-1 complex and HIF-1α-aOH, is significantly higher in DMOG treated cells than in untreated cells in normoxia. The major component of the sum in DMOG treated cells is HIF-1α, whereas its hydroxylated form predominates in normoxic cells. As expected, we observe that HIF-1β is stable in normoxia. However, in the presence of available HIF-1α, HIF-1β is consumed for HIF-1 complex formation.

FIGURE 4

Figure 4. Model simulation in NK cells without DMOG (N): (A) with one inhibitor and (B) with two inhibitors. Colors indicate the following experimental settings: (blue) untreated cells; (red) cells primed with IL-15; (green) cells primed with IL-15 and treated with STAT3 inhibitor S3I-201; (yellow) cells primed with IL-15 and treated with mTOR inhibitor rapamycin (Rapa); (magenta) cells primed with IL-15 and treated with NF-κB inhibitor (NF-κBi); (cyan) cells primed with IL-15 and treated with STAT3 inhibitor and rapamycin; (gray) cells primed with IL-15 and treated with NF-κB inhibitor and rapamycin; (orange) cells primed with IL-15 and treated with NF-κB and STAT3 inhibitors.

FIGURE 5

Figure 5. Model simulation in chemical hypoxia (20 μM DMOG): (A) with one inhibitor and (B) with two inhibitors. Colors indicate the following experimental settings: (blue) untreated cells; (red) cells primed with IL-15; (green) cells primed with IL-15 and treated with STAT3 inhibitor S3I-201; (yellow) cells primed with IL-15 and treated with mTOR inhibitor rapamycin (Rapa); (magenta) cells primed with IL-15 and treated with NF-κB inhibitor (NF-κBi); (cyan) cells primed with IL-15 and treated with STAT3 inhibitor and rapamycin; (gray) cells primed with IL-15 and treated with NF-κB inhibitor and rapamycin; (orange) cells primed with IL-15 and treated with NF-κB and STAT3 inhibitors.

Consistent with the model assumptions, we observe mTOR inhibition by rapamycin. Moreover, the stabilization of HIF-1α in DMOG treated cells and the subsequent formation of HIF-1 results in a negative feedback on mTOR (Figure 5). NF-κB shows higher activity in DMOG treated cells compared to untreated cells and our simulations predict an essential role for NF-κB as a regulator of HIF-1α-mRNA and protein in DMOG treated cells. In contrast, the IL-15-induced STAT3 activity is higher in cells without DMOG and inhibition of STAT3 results in an important reduction of HIF-1α-mRNA and protein levels. Combined inhibition of both transcription factors abolishes HIF-1α enrichment in both, DMOG treated and untreated cells.

3.4. Regulators of HIF-1α Enrichment

Parameter values used for data fitting (Table 1) in Figure 3A indicate that the external regulation rates of IL-15, mTOR, and STAT3 are small or negligible in the considered cell cultures. Nevertheless, such parameters could change from donor to donor, in particular if affected by inflammatory conditions or cancer (5, 64).

To investigate the influence of external regulators on the behavior of the total HIF-1α stabilization, we systematically perturbed the constant activation rate of different proteins in the network. First we studied the effect of external regulation of IL-15 on the stabilization of total HIF-1α in normoxia in the presence or absence of DMOG. For this we simulated ideal experiments in which the cells are exposed to continuous stimulation. We ran computer simulations varying the IL-15 external regulation parameter a₁ in the interval [0, 10] nM h⁻¹ and plotted the solution of total HIF-1α over time. All other parameter values as well as initial conditions were fixed as indicated in Materials and Methods. Figure 6 shows the results for IL-15 (left), and for total HIF-1α in the absence (middle) or presence of DMOG (right), with dark blue curves corresponding to the lowest value (a₁ = 0 nM h⁻¹) and red curves to the highest value (a₁ = 10 nM h⁻¹). Simulations confirm the synergistic effect of IL-15 and DMOG treatment on HIF-1α. The stronger the continuous IL-15 stimulus, the higher are the total HIF-1α levels. Compared with untreated cells, HIF-1α levels reach two to three times higher steady state² values in the presence of DMOG.

FIGURE 6

Figure 6. Effects of external IL-15 regulation on total HIF-1α in NK cells cultivated in the absence (N) or presence of DMOG. External IL-15 regulation rate (a₁) is varied in [0, 10] nM h⁻¹ with regular steps. Other parameters and initial values are fixed as in Tables 1, 2. Curves with same color correspond to the same parameter value and follow the jet color map in MATLAB^®, with dark blue corresponding to the lowest value (a₁ = 0 nM h⁻¹) and red to the highest value (a₁ = 10 nM h⁻¹). Left: IL-15 dynamics is equal in both, NK cells cultivated without and with DMOG as there is no feedback on IL-15 in the model (cf. Equation 2); middle: Total HIF-1α in untreated NK cells; right: Total HIF-1α in NK cells cultivated with DMOG.

Similarly, we investigated the dependence of total HIF-1α accumulation on the external regulation of mTOR and STAT3. We considered cells with or without DMOG and proceeded as above by varying the parameters a₃ (for mTOR) and a₈ (for STAT3) in the interval [0, 10] nM h⁻¹. For investigations on the steady state of the systems (t = 100 h), it makes no difference if cells are initially stimulated with IL-15 or not, as the effect of initial stimulation has vanished at the steady state (cf. Figures 4, 5). The results are shown in Figure 7, where dark blue curves correspond to the lowest value (a_j = 0 nM h⁻¹, j = 3, 8) and red curves to the highest value (a_j = 10 nM h⁻¹, j = 3, 8). Our computer simulations confirm that higher HIF-1α concentration in DMOG treated cells induces a negative feedback and downregulates mTOR (Figure 7A). As Figure 7B suggests, increasing STAT3 external regulation leads to higher HIF-1α levels, which are amplified by DMOG treatment, although STAT3 levels in DMOG treated cells are slightly lower than in cells without DMOG.

FIGURE 7

Figure 7. (A) Effects of external mTOR regulation on total HIF-1α in NK cells without (N) or with DMOG. External mTOR activation rate (a₃) is varied in [0, 10] nM h⁻¹ with regular steps. (B) Effects of external STAT3 regulation on total HIF-1α. External STAT3 activation rate (a₈) is varied in [0, 10] nM h⁻¹ with regular steps. All other parameters and initial values are fixed as in Tables 1, 2. Curves with same color correspond to the same parameter value and follow the jet color map in MATLAB^®, with dark blue corresponding to the lowest value (a_j = 0 nM h⁻¹, j = 3, 8) and red to the highest value (a_j = 10 nM h⁻¹, j = 3, 8).

We further investigated the dependence of total HIF-1α enrichment on two signals at the same time. We started by varying the external activation rate of mTOR (a₃) and STAT3 (a₈) in the interval [0, 10] nM h⁻¹, ran simulations up to t = 100 h and obtained numerical solutions of the mathematical model at the steady state. Figure 8 shows the results for total HIF-1α, STAT3, and mTOR in NK cells cultivated with or without DMOG. The same figure shows also the effect of simultaneous changes in the external activation rate of NF-κB (a₇ in the interval [0, 10] nM h⁻¹) and STAT3 (a₈ in the interval [0, 10] nM h⁻¹). Again, we observe a central role of STAT3 as regulator of HIF-1α enrichment, especially in synergy with NF-κB. An increase in STAT3 external regulation rate combined with an increase in the external regulation rate of NF-κB leads to a higher amplification of total HIF-1α compared to STAT3 combined with mTOR. Plots in Figure 8 also reflect the negative feedback of induced HIF-1 on mTOR: (i) in general, NK cells treated with DMOG have lower levels of activated mTOR than untreated cells and (ii) higher concentration of activated STAT3 (due to increasing a₈ rate) induces HIF-1α, resulting in higher levels of HIF-1 complex, which in turn is known to inhibit mTOR. Finally, Figure 8 stresses the role of HIF-1 as activator of NF-κB. In normoxic cells, where HIF-1 levels are low, NF-κB activity is low, despite increasing of external regulation. In contrast, in DMOG treated cells, HIF-1 accumulation leads to upregulation of NF-κB. This is in accordance with data obtained in neutrophils demonstrating that NF-κB is an important downstream effector of the HIF-1α-dependent response (37). Figure 9 shows HIF-1α steady states in dependence on the external regulation rate of IL-15 (a₁ varying in [0, 10] nM h⁻¹) and activation of STAT3 (a₈ varying in the interval [0, 10] nM h⁻¹). The results confirm the synergistic effect of IL-15, STAT3 and DMOG in increasing HIF-1α levels.

FIGURE 8

Figure 8. Effects of external regulation of mTOR, NF-κB, and STAT3 in NK cells without (N) or with DMOG. External activation rates for mTOR (a₃) and STAT3 (a₈) are varied in the interval [0, 10] nM h⁻¹ and for NF-κB (a₇) in the interval [10, 20] nM h⁻¹. Other parameters and initial values are fixed as in Tables 1, 2. Steady states (100 h) of the model solutions are computed for total HIF-1α (first row), STAT3 (second row), and mTOR or NF-κB (third row).

FIGURE 9

Figure 9. Effects of external regulation of IL-15 (a₁) and STAT3 (a₈) in NK cells without DMOG (N, left) and DMOG treated cells (right). External regulation rates are varied in the intervals [0, 5] nM h⁻¹ for IL-15 and [0, 10] nM h⁻¹ for STAT3. Other parameters and initial values are fixed as in Tables 1, 2. Steady states (100 h) of the model solutions are computed for total HIF-1α (first row) and STAT3 (second row).

Besides the above deterministic perturbations, we tested the network robustness with a stochastic approach. Robustness allows a system to maintain its function, regardless of external and internal perturbations (65). We perturbed specific elements of the system 25% around their originally estimated value (parameter values are otherwise fixed as in Tables 1, 2). We applied stochastic perturbations using the Monte Carlo method (see section 2.8) and computed how external (in IL-15) and internal (in mTOR and AKT) changes affect the steady state of total HIF-1α in IL-15 stimulated cells treated with DMOG. Figure 10 shows the histograms for IL-15 (A), AKT (B), and mTOR (C) (green) and total HIF-1α (red). Our results show that the model is robust to internal and external stochastic perturbations, indicating that variations of ± 25% in IL-15, in the AKT activation rate a₂ ± 25% or in the mTOR activation rate a₃ ± 25% result in minimal variations (<10%) in the steady state of total HIF-1α.

FIGURE 10

Figure 10. Stochastic extrinsic and intrinsic perturbations of the HIF-1α regulatory network in IL-15 stimulated cells treated with DMOG. Histograms accounting for total HIF-1α relative steady state variations after stochastically varying (A) IL-15, (B) the AKT basal activation rate (a₂), and (C) the mTOR basal activation rate (a₃). The model is robust to internal and external stochastic perturbations, indicating that variation of ± 25% in the above elements result in minimal variations (<10%) in the steady state of total HIF-1α.

With the help of numerical simulations we tested how HIF-1α stabilization is affected by increasing concentration of DMOG. Assuming that NK cells are treated with DMOG and stimulated with IL-15 at time t = 0 h, we changed the DMOG concentration from 0 to 100%, with 100% corresponding to 20 μM. Figure 11A shows the evolution of HIF-1α in time, with HIF-1α stabilization depending on the DMOG dosage. We computed the fold change of HIF-1α stabilization at the equilibrium (t = 100 h) and compared control cells (untreated) with cells treated with different DMOG concentrations (Figure 11B). The relation between HIF-1α stabilization and DMOG dosage is non-linear and doubling the DMOG dose does not lead to twice as high HIF-1α levels. Our results suggest an exponential trend in the relation between PHD/FIH inhibitor DMOG and HIF-1α stabilization, which is in accordance with what was previously observed for HEK cells (19).

FIGURE 11

Figure 11. Stabilization of HIF-1α in dependence on DMOG concentration after IL-15 stimulation. (A) Evolution of HIF-1α in time, depending on the DMOG dosage with dark blue corresponding to no DMOG (N) and red corresponding to 20 μM DMOG. In these simulations cells are initially stimulated with IL-15. (B) Fold change of HIF-1α stabilization at the equilibrium (t = 100 h). Simulation of NK cell treatment with increasing concentrations of DMOG. Parameter values for these simulations are chosen as in Tables 1, 2, with 100% corresponding to 20 μM DMOG.

3.5. Which Timing for Cell Treatment?

Model (2)–(11) can be used for a number of in silico experiments to test the validity of biological hypotheses or predict the outcome of laboratory tests. In this study we were particularly interested in the synergy of IL-15-stimulation and DMOG treatment in the stabilization of HIF-1α, already observed in the results described above. In all previous simulations, normoxic NK cells were stimulated at the beginning of the observation with IL-15 in the presence or absence of DMOG. To understand how the timing of the treatments affects HIF-1α stabilization in NK cells, we also simulated different possibilities for the timing of cell treatment combining chemical hypoxia and stimulation with IL-15 (Figure 12).

FIGURE 12

Figure 12. (A) HIF-1α evolution in time and (B) fold change at t = 12 h, normalized with respect to HIF-1α levels in untreated cells. Colors correspond to the following experimental settings: (gray) untreated cells (N) for 30 h; (dark blue) IL-15 stimulation at time t = 0 h; (magenta) cells in DMOG for 30 h; (green) cells in DMOG for 30 h with IL-15 stimulation at t = 0 h; (yellow) IL-15 stimulation at t = 6 h; (orange) DMOG treatment at t = 6 h; (red) IL-15 stimulation and DMOG treatment at t = 6 h; (light blue) DMOG treatment at t = 0 h and IL-15 stimulation at t = 6 h; (black) IL-15 stimulation at t = 0 h and DMOG treatment at t = 6 h.

We compared the HIF-1α dynamics for the following in silico experiments: (gray) untreated NK cells (N) cultivated for 30 h; (dark blue) IL-15 stimulation at t = 0 h; (magenta) DMOG treatment for 30 h; (green) DMOG treatment for 30 h, with IL-15 stimulation at t = 0 h; (yellow) IL-15 stimulation at t = 6 h; (orange) DMOG treatment starting at t = 6 h; (red) IL-15 stimulation at t = 6 h and DMOG treatment starting at t = 6 h; (light blue) DMOG treatment for 30 h with IL-15 stimulation at t = 6 h; (black) IL-15 stimulation at t = 0 h and DMOG treatment starting at t = 6 h.

Figure 12A shows the time evolution of HIF-1α stabilization over 30 h. We observe the impulses at t = 6 h due to changes in the experimental conditions. On the long term, the effect of IL-15 stimulation vanishes and HIF-1α levels converge to those reached in unstimulated cells. Figure 12B shows the fold change of HIF-1α at t = 12 h. Values are normalized with respect to HIF-1α in untreated cells (N, gray bar). We observe that on the short time scale the timing of treatments importantly affects HIF-1α stabilization. In particular, treating the cells first with DMOG or first stimulating them with IL-15 is not equivalent (compare the black bar and the light blue bar). The highest HIF-1α levels after 12 h are reached when cells are first stimulated with IL-15 at t = 0 and treated with DMOG at t = 6 h (black bar). Cultivating cells in normoxia and treating them with IL-15 and DMOG at time t = 6 h (red) yields lower HIF-1α values than 12 h cultivation in the presence of DMOG after initial (t = 0 h) IL-15 stimulation (green).

4. Discussion

Being an essential mediator of cellular adaptation to hypoxia (66, 67), HIF-1α plays a critical role as regulator of inflammation and immune system response (36, 68). The understanding of its regulation is crucial in immunology.

While HIF-1α hydroxylation and degradation pathways have been recently described using mathematical methods (19, 20), less is known concerning the mechanistic description of processes regulating the levels of HIF-1α mRNA and protein (10). In this work we have presented a combined approach of experimental and mathematical analysis to understand HIF-1α regulation in human NK cells, in particular simulating hypoxic (DMOG) and inflammatory (IL-15) conditions. To the best of our knowledge, there is no previous interdisciplinary approach describing the interplay of hypoxia and IL-15 stimulation, and their effects on HIF-1α dynamics in immune cells. The proposed mathematical model (2)–(11) and the estimated parameter values (Tables 1, 2) explain collected time series for HIF-1α, also catching the dynamics of other regulatory proteins (Figure 3). Our simulation results and in silico experiments highlight the synergy of IL-15 and hypoxia in HIF-1α stabilization, suggesting an important role for STAT3 and NF-κB as regulators of IL-15 induced HIF-1α enrichment in peripheral NK cells.

The mathematical model proposed in this work aimed at the qualitative mechanistic description of IL-15 induced biochemical processes regulating HIF-1α stabilization in NK cells. We made use of collected time series (HIF-1α, AKT, and STAT3) for quantitative investigation, data fitting and model predictions. A limitation of our results is that model predictions for quantities lacking experimental information (e.g., NF-κB, mTOR, and HIF-1α-mRNA) can be made only on a relative scale (43). While the calibration (Figure 3A) and validation results (Figure 3B) for HIF-1α are overall very satisfactory, the quantitative match for STAT3 in cells treated with DMOG and IL-15 (Figure 3B) could be improved. This might be achieved by refining the fit for STAT3 time series in IL-15-stimulated NK cells treated with DMOG and rapamycin (Figure 3A). In this study, we performed an all-at-once parameter estimation, applying a direct multiple shooting approach and a gradient-based (generalized Gauss-Newton) method (cf. section 2.7). The method is known to perform well and converge fast [cf. (45)] but, being a local optimization method, it might get stuck in a local optimal minimum. A possible method to overcome local minima is to perform many independent optimization runs starting from randomly selected starting points (69). Alternatively, one could adopt global optimization methods, which however can be computationally very costly (45, 69).

Concerning the statistical significance of the parameter estimates (Table 2) we have adopted here a first order approximation of non-linear confidence regions. Parameter estimation and identifiability could be further refined and investigated, e.g., performing a second-order analysis of the non-linear confidence regions [cf. (70)] or exploiting the profile likelihood, as suggested by Raue et al. (69).

Our model captures essential features of HIF-1α regulation, making a number of simplifying assumptions. Several model extensions and refinements could be proposed. For example, we could include further steps in the degradation pathway of HIF-1α, as proposed by others (10, 19). Moreover, the dynamics of IL-15 is simply given by constant production and degradation rates [as it has previously been assumed by other authors, e.g., for IL-21 dynamics (71)], and sensitivity analysis indicates that the AKT dynamics is approximatively linear in the considered experimental setting (section 3.3). The reaction cascade downstream of IL-15 involves several components, including the IL-15Rβγ-subunits (4), which are known to be constitutively expressed on NK cells (72) but were neglected in the proposed mathematical model. Further, the IL-15-induced activation of AKT, NF-κB and STAT3 is modeled by means of linear terms. One possible model extension would include non-linear terms (Michaelis-Menten or higher order Hill functions) for the activation of IL-15 regulated proteins. Factors connected to IL-15 stimulations, such as IL-15 receptor binding and trafficking or other IL-15 induced signaling cascades (JAK/STAT5, Ras-Raf-MEK), might affect NK cell response to this cytokine and could also be taken into consideration. Further, the relation between mTOR and the HIF-1 complex could be investigated in detail. We have assumed that hypoxia downregulates mTOR (33) by means of a negative feedback of HIF-1 on mTOR. Nonetheless, the regulatory mechanism of mTOR is far more complicated, involving REDD1 and the Tsc1/Tsc2 complex (73). Our experimental data and modeling results show that HIF-1α accumulation in cells stimulated with IL-15 and treated with DMOG correlates with reduction of STAT3 activity. Our modeling approach suggests (cf. section S2) that the known negative feedback of HIF-1 on mTOR (33) is amplified by a direct inhibitory effect of DMOG on IL-15-induced STAT3 activation. This means that the observed STAT3 inhibition could not only be due to chemical hypoxia, stabilizing HIF-1α, but also due to additional effects of DMOG on NK cells. To further explore the role of DMOG on IL-15-induced STAT3 activation in NK cells, the experiments proposed in this study could be performed in cells cultivated in hypoxia (1% O₂) instead of chemical hypoxia.

Spatial effects could also be taken into account for model refinement. In contrast to previous studies (19, 74), in our modeling approach we did not make any distinction between proteins in the cell cytoplasm and the nucleus, but simply consider total cellular concentrations. In general, increased model complexity necessarily calls for more detailed experimental data in order to achieve adequate model calibration and trustworthy predictions.

We hypothesize that the proposed regulatory network is appropriate for describing HIF-1α regulation not only in NK cells but also in other types of immune cells. Moreover, the model (2)–(11) can be applied to refine and extend mathematical models in which HIF-1α dynamics is involved, e.g., models of cell cycle regulation (75) or cell proliferation (76). When studying the effects of biochemical signaling at the cellular level, it might be convenient to adopt simpler regulatory networks than those proposed here [see for example (77) for a model for proliferation of IL-15 stimulated NK cells]. To this extent model reduction could be performed by mean of biological assumptions, e.g., via quasi-steady state approximations. Further, (global) sensitivity analysis results could be used to rank the relative influence of the model parameters on the model output, and could suggest how to simplify the regulatory network, identifying parameters that minimally impact model outputs [cf. (78, 79)].

Being involved in cytokine expression, myeloid cell migration and effector functions, HIF-1α regulates both innate and adaptive immunity (80). Understanding the molecular mechanisms involved in the regulation of the HIF-1α pathway, in particular in immune cells, is of central importance to the immune cell function and could be a promising strategy in the design of treatments for human inflammatory diseases and cancer. Our results indicate that NF-κB and STAT3 are important regulators of HIF-1α enrichment in IL-15 stimulated NK cells. It is tempting to speculate that a secondary effect of pharmacological STAT3 inhibition in cancer therapy may consist in a reduction of IL-15 dependent HIF-1α enrichment in NK cells, which may be expected to improve NK cell anti-tumor activity (12, 13).

Data Availability Statement

All datasets generated for this study are included in the manuscript/Supplementary Files.

Ethics Statement

The study was reviewed and approved by the Medical Ethics Commission II of the Medical Faculty Mannheim, Heidelberg University (2014-500N-MA).

Author Contributions

MB and HL: conceptualization. AC, MB, EK, and HL: methodology. SV: NK cell characterization. AB, MB, and AF: software. AF: robustness analysis. AC and MB: investigation and data curation. MB and AC: writing. AC, MB, SV, MT, and HL: review and editing (original manuscript). MB, AC, AB, AF, EK, and H-GB: review and editing (revision). All authors approved the final version of the manuscript.

Funding

This research work was supported by funds to H-GB, EK, HL, and MT by the Klaus Tschira Foundation, Germany. MB was supported by the European Social Fund and by the Ministry of Science, Research and Arts Baden-Württemberg.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Jutta Schulte and Bianca S. Himmelhan (Department of Anesthesiology and Surgical Intensive Care Medicine, University Medical Center Mannheim, Medical Faculty Mannheim, Heidelberg University) for technical support. Further we acknowledge financial support by Deutsche Forschungsgemeinschaft within the funding program Open Access Publishing, by the Baden-Württemberg Ministry of Science, Research and the Arts and by Ruprecht-Karls-Universität Heidelberg.

The authors would like to thank the referees for their valuable critical comments which helped to improve the quality of this manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2019.02401/full#supplementary-material

Footnotes

1. ^The dynamics of IL-15 is trivial and neither affected by oxygen saturation nor other inhibitors (in the model assumptions there is no feedback on IL-15, cf. Equation 2), hence it is equivalent in all experimental conditions and omitted in Figures 4, 5. The same holds for the AKT dynamics in view of the results of the sensitivity analysis.

2. ^A complete analytical study of the steady states of system (2)–(11) and their stability was beyond the scope of this manuscript. In general, the non-linear system (2)–(11) might have, in certain parameter ranges, multiple biologically relevant steady states. However, the parameter values used for the numerical simulations in this work guarantee convergence to a unique non-negative steady state which is shifted to higher/lower values when different stimuli are considered (cf. simulations shown in Figures 4–10).

References

1. Chan C, Smyth M, Martinet L. Molecular mechanisms of natural killer cell activation in response to cellular stress. Cell Death Differ. (2014) 21:5. doi: 10.1038/cdd.2013.26

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Long EO, Sik Kim H, Liu D, Peterson ME, Rajagopalan S. Controlling natural killer cell responses: integration of signals for activation and inhibition. Annu Rev Immunol. (2013) 31:227–58. doi: 10.1146/annurev-immunol-020711-075005

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Mishra A, Sullivan L, Caligiuri MA. Molecular pathways: interleukin-15 signaling in health and in cancer. Clin Cancer Res. (2014) 20:2044–50. doi: 10.1158/1078-0432.CCR-12-3603

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Van den Bergh JM, Van Tendeloo VF, Smits EL. Interleukin-15: new kid on the block for antitumor combination therapy. Cytokine Growth Factor Rev. (2015) 26:15–24. doi: 10.1016/j.cytogfr.2014.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Rautela J, Huntington ND. IL-15 signaling in NK cell cancer immunotherapy. Curr Opin Immunol. (2017) 44:1–6. doi: 10.1016/j.coi.2016.10.004

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Mandal A, Viswanathan C. Natural killer cells: in health and disease. Hematol Oncol Stem Cell Ther. (2015) 8:47–55. doi: 10.1016/j.hemonc.2014.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Fogler WE, Volker K, McCormick KL, Watanabe M, Ortaldo JR, Wiltrout RH. NK cell infiltration into lung, liver, and subcutaneous B16 melanoma is mediated by VCAM-1/VLA-4 interaction. J Immunol. (1996) 156:4707–14.

PubMed Abstract | Google Scholar

8. Glas R, Franksson L, Une C, Eloranta ML, Öhlén C, Örn A, et al. Recruitment and activation of natural killer (NK) cells in vivo determined by the target cell phenotype: an adaptive component of NK cell–mediated responses. J Exp Med. (2000) 191:129–38. doi: 10.1084/jem.191.1.129

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Semenza GL. Hypoxia-inducible factors in physiology and medicine. Cell. (2012) 148:399–408. doi: 10.1016/j.cell.2012.01.021

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Fábián Z, Taylor CT, Nguyen LK. Understanding complexity in the HIF signaling pathway using systems biology and mathematical modeling. J Mol Med. (2016) 94:377–390. doi: 10.1007/s00109-016-1383-6

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Palazón A, Goldrath AW, Nizet V, Johnson RS. HIF transcription factors, inflammation, and immunity. Immunity. (2014) 41:518–28. doi: 10.1016/j.immuni.2014.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Krzywinska E, Kantari-Mimoun C, Kerdiles Y, Sobecki M, Isagawa T, Gotthardt D, et al. Loss of HIF-1α in natural killer cells inhibits tumour growth by stimulating non-productive angiogenesis. Nat Commun. (2017) 8:1597. doi: 10.1038/s41467-017-01599-w

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Ni J, Bühler L, Stojanovic A, Arnold A, Sexl V, Cerwenka A. Inhibition of the HIF-1α-mediated checkpoint refuels NK activation in cancer [abstract]. In: Proceedings of the American Association for Cancer Research Annual Meeting 2018; 2018 Apr 14-18; Chicago, IL. Philadelphia (PA): AACR. Cancer Res. (2018) 78:Abstract nr 4743. doi: 10.1158/1538-7445.AM2018-4743

CrossRef Full Text | Google Scholar

14. Kohn KW, Riss J, Aprelikova O, Weinstein JN, Pommier Y, Barrett JC. Properties of switch-like bioregulatory networks studied by simulation of the hypoxia response control system. Mol Biol Cell. (2004) 15:3042–52. doi: 10.1091/mbc.e03-12-0897

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Kooner P, Maini PK, Cavaghan D. Mathematical modeling of the HIF-1 mediated hypoxic response in tumours. In: Proceedings of the 2005 International Symposium on Mathematical and Computational Biology. Rio de Janeiro: E-papers (2006). p. 281–315.

Google Scholar

16. Qutub AA, Popel AS. A computational model of intracellular oxygen sensing by hypoxia-inducible factor HIF-1α. J Cell Sci. (2006) 119:3467–80. doi: 10.1242/jcs.03087

CrossRef Full Text | Google Scholar

17. Dayan F, Monticelli M, Pouysségur J, Pécou E. Gene regulation in response to graded hypoxia: the non-redundant roles of the oxygen sensors PHD and FIH in the HIF pathway. J Theor Biol. (2009) 259:304–16. doi: 10.1016/j.jtbi.2009.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Schmierer B, Novák B, Schofield CJ. Hypoxia-dependent sequestration of an oxygen sensor by a widespread structural motif can shape the hypoxic response–a predictive kinetic model. BMC Syst Biol. (2010) 4:139. doi: 10.1186/1752-0509-4-139

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Nguyen LK, Cavadas MA, Scholz CC, Fitzpatrick SF, Bruning U, Cummins EP, et al. A dynamic model of the Hypoxia-Inducible Factor 1α (HIF-1α) network. J Cell Sci. (2013) 126:1454–63. doi: 10.1242/jcs.119974

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Cavadas MA, Nguyen LK, Cheong A. Hypoxia-inducible factor (HIF) network: insights from mathematical models. Cell Comm Signal. (2013) 11:42. doi: 10.1186/1478-811X-11-42

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Velásquez SY, Killian D, Schulte J, Sticht C, Thiel M, Lindner HA. Short-term hypoxia synergizes with interleukin 15 priming in driving glycolytic gene transcription and supports human natural killer cell activities. J Biol Chem. (2016) 291:12960–77. doi: 10.1074/jbc.M116.721753

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Brauer F, Nohel JA. The Qualitative Theory of Ordinary Differential Equations: An Introduction. New York, NY: Dover Publication (1989).

Google Scholar

23. Müller J, Kuttler C. Methods and Models in Mathematical Biology. Berlin; Heidelberg: Springer (2015).

Google Scholar

24. Marçais A, Cherfils-Vicini J, Viant C, Degouve S, Viel S, Fenis A, et al. The metabolic checkpoint kinase mTOR is essential for IL-15 signaling during the development and activation of NK cells. Nat Immunol. (2014) 15:749. doi: 10.1038/ni.2936

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Nandagopal N, Alaa KA, Amandeep KK, Lee SH. The critical role of IL-15-PI3K-mTOR pathway in natural killer cell effector functions. Front Immunol. (2014) 5:187. doi: 10.3389/fimmu.2014.00187

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Wagner JA, Rosario M, Romee R, Berrien-Elliott MM, Schneider SE, Leong JW, et al. CD56bright NK cells exhibit potent antitumor responses following IL-15 priming. J Clin Invest. (2017) 127:4042–58. doi: 10.1172/JCI90387

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Fehniger TA, Caligiuri MA. Interleukin 15: biology and relevance to human disease. Blood. (2001) 97:14–32. doi: 10.1182/blood.V97.1.14

PubMed Abstract | CrossRef Full Text | Google Scholar

28. McDonald PP, Russo MP, Ferrini S, Cassatella MA. Interleukin-15 (IL-15) induces NF-κB activation and IL-8 production in human neutrophils. Blood. (1998) 92:4828–35.

PubMed Abstract | Google Scholar

29. Budagian V, Bulanova E, Paus R, Bulfone-Paus S. IL-15/IL-15 receptor biology: a guided tour through an expanding universe. Cytokine Growth Factor Rev. (2006) 17:259–80. doi: 10.1016/j.cytogfr.2006.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Kelly E, Won A, Refaeli Y, Parijs V. IL-2 and related cytokines can promote T cell survival by activating AKT. J Immunol. (2002) 168:597–603. doi: 10.4049/jimmunol.168.2.597

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Xu Q, Briggs J, Park S, Niu G, Kortylewski M, Zhang S, et al. Targeting STAT3 blocks both HIF-1 and VEGF expression induced by multiple oncogenic growth signaling pathways. Oncogene. (2005) 24:5552. doi: 10.1038/sj.onc.1208719

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Pouysségur J, Dayan F, Mazure NM. Hypoxia signalling in cancer and approaches to enforce tumour regression. Nature. (2006) 441:437–43. doi: 10.1038/nature04871

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Brugarolas J, Lei K, Hurley RL, Manning BD, Reiling JH, Hafen E, et al. Regulation of mTOR function in response to hypoxia by REDD1 and the TSC1/TSC2 tumor suppressor complex. Genes Dev. (2004) 18:2893–904. doi: 10.1101/gad.1256804

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Weichhart T, Hengstschläger M, Linke M. Regulation of innate immune cell function by mTOR. Nat Rev Immunol. (2015) 15:599–614. doi: 10.1038/nri3901

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Dan HC, Cooper MJ, Cogswell PC, Duncan JA, Ting JP, Baldwin AS. Akt-dependent regulation of NF-κB is controlled by mTOR and raptor in association with IKK. Genes Dev. (2008) 22:1490–500. doi: 10.1101/gad.1662308

CrossRef Full Text | Google Scholar

36. D'Ignazio L, Bandarra D, Rocha S. NF-κB and HIF crosstalk in immune responses. FEBS J. (2016) 283:413–24. doi: 10.1111/febs.13578

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Walmsley SR, Farahi N, Peyssonnaux C, Johnson RS, Cramer T, Sobolewski A, et al. Hypoxia-induced neutrophil survival is mediated by HIF-1α dependent NF-κB activity. J Exp Med. (2005) 201:105–15. doi: 10.1084/jem.20040624

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Lancaster DE, McDonough MA, Schofield CJ. Factor inhibiting hypoxia-inducible factor (FIH) and other asparaginyl hydroxylases. Biochem Soc Trans. (2004) 32:943–5. doi: 10.1042/BST0320943

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Shampine LF, Reichelt MW. The MATLAB ODE suite. SIAM J Scie Comp. (1997) 18:1–22. doi: 10.1137/S1064827594276424

CrossRef Full Text | Google Scholar

40. Bauer I, Bock HG, Schlöder JP. DAESOL–A BDF-Code for the Numerical Solution of Differential Algebraic Equations (Version 3.0.1). Heidelberg: IWR Heidelberg University (1999).

41. Bauer I, Bock HG, Körkel S, Schlöder JP. Numerical methods for initial value problems and derivative generation for DAE models with application to optimum experimental design of chemical processes. In: Keil F, Mackens W, VoßH, Werther J, editors. Scientific Computing in Chemical Engineering II. Berlin; Heidelberg: Springer (1999). p. 282–9.

Google Scholar

42. Körkel S. Numerische Methoden für Optimale Versuchsplanungsprobleme bei nichtlinearen DAE-Modellen (Dissertation), Universität Heidelberg, Heidelberg, Germany (2002). Available online at: http://www.koerkel.de

43. Raue A, Schilling M, Bachmann J, Matteson A, Schelke M, Kaschek D, et al. Lessons learned from quantitative dynamical modeling in systems biology. PLoS ONE. (2013) 8:e74335. doi: 10.1371/journal.pone.0074335

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Fiedler A, Raeth S, Theis FJ, Hausser A, Hasenauer J. Tailored parameter optimization methods for ordinary differential equation models with steady-state constraints. BMC Syst Biol. (2016) 10:80. doi: 10.1186/s12918-016-0319-7

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Ashyraliyev M, Fomekong-Nanfack Y, Kaandorp JA, Blom JG. Systems biology: parameter estimation for biochemical models. FEBS J. (2009) 276:886–902. doi: 10.1111/j.1742-4658.2008.06844.x

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Bard Y. Nonlinear Parameter Estimation. New York, NY: Academic Press (1974).

Google Scholar

47. Bates DM, Watts DG. Nonlinear Regression Analysis and Its Applications. Wiley series in probability and mathematical statistics: Applied probability and statistics. Chichester: Wiley (1988).

Google Scholar

48. López C D, Barz T, Körkel S, Wozny G. Nonlinear ill-posed problem analysis in model-based parameter estimation and experimental design. Comput Chem Eng. (2015) 77:24–42. doi: 10.1016/j.compchemeng.2015.03.002

CrossRef Full Text | Google Scholar

49. Bock HG, Körkel S, Kostina E, Schlöder JP. Robustness aspects in parameter estimation, optimal design of experiments and optimal control. In: Jäger W, Rannacher R, Warnatz J, editors. Reactive Flows, Diffusion and Transport. From Experiments via Mathematical Modeling to Numerical Simulation and Optimization: Final Report of SFB (Collaborative Research Center) 359. Heidelberg: Springer (2007). p. 117–46.

Google Scholar

50. Bock HG, Kostina E, Schlöder JP. Direct multiple shooting and generalized Gauss-Newton method for parameter estimation problems in ODE models. In: Carraro T, Geiger M, Körkel S, Rannacher R, editors. Multiple Shooting and Time Domain Decomposition Methods. Cham: Springer International Publishing (2015). p. 1–34.

Google Scholar

51. Bock HG. Randwertproblemmethoden zur Parameteridentifizierung in Systemen Nichtlinearer Differentialgleichungen. 183. Bonn: Der Math.-Naturwiss. Fakultät der Universität Bonn (1987).

52. Bock HG, Kostina E, Schlöder JP. Numerical Methods for Parameter Estimation in Nonlinear Differential Algebraic Equations. GAMM Mitteilungen. 2007;30/2:376–408.

Google Scholar

53. Bock HG. Numerical treatment of inverse problems in chemical reaction kinetics. In: Ebert KH, Deuflhard P, Jäger W, editors. Modelling of Chemical Reaction Systems. Vol. 18 of Springer Series in Chemical Physics. Heidelberg: Springer (1981). p. 102–25. Available online at: https://link.springer.com/chapter/10.1007/978-3-642-68220-9_8

Google Scholar

54. Schlöder JP. Numerische Methoden zur Behandlung Hochdimensionaler Aufgaben der Parameteridentifizierung. Vol. 187 of Bonner Mathematische Schriften. Bonn: Universität Bonn (1988).

55. Allgower EL, Georg K. Numerical Continuation Methods: An Introduction. Springer Science & Business Media (1990).

Google Scholar

56. Olufsen MS, Ottesen JT. A practical approach to parameter estimation applied to model predicting heart rate regulation. J Math Biol. (2013) 67:39–68. doi: 10.1007/s00285-012-0535-8

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Pope SR. Parameter identification in lumped compartment cardiorespiratory models (PhD thesis). Applied Mathematics, NC State University, Raleigh, NC, United States (2009).

Google Scholar

58. Semenza GL. Targeting HIF-1 for cancer therapy. Nat Rev Cancer. (2003) 3:721. doi: 10.1038/nrc1187

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Palazón A, Aragonés J, Morales-Kastresana A, de Landázuri MO, Melero I. Molecular pathways: hypoxia response in immune cells fighting or promoting cancer. Clin Cancer Res. (2012) 18:1207–13. doi: 10.1158/1078-0432.CCR-11-1591

CrossRef Full Text | Google Scholar

60. Dang EV, Barbi J, Yang HY, Jinasena D, Yu H, Zheng Y, et al. Control of T H 17/T reg balance by hypoxia-inducible factor 1. Cell. (2011) 146:772–84. doi: 10.1016/j.cell.2011.07.033

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Meng X, Grötsch B, Luo Y, Knaup KX, Wiesener MS, Chen XX, et al. Hypoxia-inducible factor-1α is a critical transcription factor for IL-10-producing B cells in autoimmune disease. Nat Commun. (2018) 9:251. doi: 10.1038/s41467-017-02683-x

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Alessi DR, Andjelkovic M, Caudwell B, Cron P, Morrice N, Cohen P, et al. Mechanism of activation of protein kinase B by insulin and IGF-1. EMBO J. (1996) 15:6541–51. doi: 10.1002/j.1460-2075.1996.tb01045.x

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Wen Z, Zhong Z, Darnell JE Jr. Maximal activation of transcription by STATl and STAT3 requires both tyrosine and serine phosphorylation. Cell. (1995) 82:241–50. doi: 10.1016/0092-8674(95)90311-9

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Yu H, Pardoll D, Jove R. STATs in cancer inflammation and immunity: a leading role for STAT3. Nat Rev Cancer. (2009) 9:798. doi: 10.1038/nrc2734

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Kitano H. Biological robustness. Nat Rev Genet. (2004) 5:826. doi: 10.1038/nrg1471

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Semenza GL. Regulation of mammalian O2 homeostasis by hypoxia-inducible factor 1. Annu Rev Cell Dev Biol. (1999) 15:551–78. doi: 10.1146/annurev.cellbio.15.1.551

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Sitkovsky M, Lukashev D. Regulation of immune cells by local-tissue oxygen tension: HIF-1α and adenosine receptors. Nat Rev Immunol. (2005) 5:712. doi: 10.1038/nri1685

CrossRef Full Text | Google Scholar

68. Nizet V, Johnson RS. Interdependence of hypoxic and innate immune responses. Nat Rev Immunol. (2009) 9:609. doi: 10.1038/nri2607

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Raue A, Kreutz C, Maiwald T, Bachmann J, Schilling M, Klingmüller U, et al. Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics. (2009) 25:1923–29. doi: 10.1093/bioinformatics/btp358

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Kostina E, Nattermann M. Second-order sensitivity analysis of parameter estimation problems. Int J Uncert Quant. (2015) 5:209–331. doi: 10.1615/Int.J.UncertaintyQuantification.2015010312

CrossRef Full Text | Google Scholar

71. Cappuccio A, Elishmereni M, Agur Z. Cancer immunotherapy by interleukin-21: potential treatment strategies evaluated in a mathematical model. Cancer Res. (2006) 66:7293–300. doi: 10.1158/0008-5472.CAN-06-0241

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Carson WE, Fehniger TA, Haldar S, Eckhert K, Lindemann MJ, Lai CF, et al. A potential role for interleukin-15 in the regulation of human natural killer cell survival. J Clin Invest. (1997) 99:937–43. doi: 10.1172/JCI119258

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Caron E, Ghosh S, Matsuoka Y, Ashton-Beaucage D, Therrien M, Lemieux S, et al. A comprehensive map of the mTOR signaling network. Mol Syst Biol. (2010) 6:453. doi: 10.1038/msb.2010.108

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Bedessem B, Stéphanou A. Role of compartmentalization on HIF-1α degradation dynamics during changing oxygen conditions: a computational approach. PLoS ONE. (2014) 9:e110495. doi: 10.1371/journal.pone.0110495

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Bedessem B, Stéphanou A. A mathematical model of HIF-1α-mediated response to hypoxia on the G1/S transition. Math Biosci. (2014) 248:31–9. doi: 10.1016/j.mbs.2013.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Zhang B, Ye H, Yang A. Mathematical modelling of interacting mechanisms for hypoxia mediated cell cycle commitment for mesenchymal stromal cells. BMC Syst Biol. (2018) 12:35. doi: 10.1186/s12918-018-0560-3

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Zhao YM, French AR. Two-compartment model of NK cell proliferation: insights from population response to IL-15 stimulation. J Immunol. (2012) 188:2981–90. doi: 10.4049/jimmunol.1102989

PubMed Abstract | CrossRef Full Text | Google Scholar

78. van Riel NA. Dynamic modelling and analysis of biochemical networks: mechanism-based models and model-based experiments. Briefings Bioinf. (2006) 7:364–74. doi: 10.1093/bib/bbl040

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Wentworth MT, Smith RC, Banks HT. Parameter selection and verification techniques based on global sensitivity analysis illustrated for an HIV model. SIAM/ASA J Uncert Quant. (2016) 4:266–97. doi: 10.1137/15M1008245

CrossRef Full Text | Google Scholar

80. Imtiyaz HZ, Simon MC. Hypoxia-inducible factors as essential regulators of inflammation. Curr Top Microbiol Immunol. (2010) 345:105–20. doi: 10.1007/82_2010_74

PubMed Abstract | CrossRef Full Text | Google Scholar