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.


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)(7)(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)(15)(16)(17)(18)(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.

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.

In vitro Treatments
Freshly isolated NK cells were maintained overnight under standard conditions and were stimulated with human recombinant IL-15 (

Western Blotting
Total cell extracts were prepared by resuspending 3×10 6 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.

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.

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 with states y, controls u and parameters p. The vector y(t) ∈ R n y indicates the "state of the system, " that is, the concentration of the considered proteins, complexes and mRNA at time t ∈ [t 0 , t f ] ⊂ R, and y 0 ∈ R n y is the initial state. The parameter vector p ∈ R 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) ∈ R 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 R 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.
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)(28)(29). All in all, we assumed that IL-15 (y 1 ) activates AKT (y 2 ), NF-κB (y 7 ), and STAT3 (y 8 ). For a general formulation we further assumed that IL-15 enters the system at constant rate a 1 and decays at rate d 1 , that is, the IL-15 dynamics is given by Activation of AKT is assumed to occur via IL-15 (30) (activation rate k 1 ), other external mediators (basal activation rate a 2 ), and also via STAT3 (maximal rate k S ) (31). We further assumed that AKT constantly decays at rate d 2 , yielding y ′ 2 (t) = a 2 + k 1 y 1 + k S y n 2 8 ξ n 2 28 + y n 2 Following the literature (32) we assumed that AKT activates mTOR (y 3 ) at rate k 2 . mTOR basal activation and decay rate are denoted by a 3 and d 3 , respectively. The inhibitory effect that hypoxia has on mTOR (33) is included in the model by means of In the left column we report the parameter value used for the numerical simulations. In the right column we indicate whether the parameter value has been fixed because of previous literature, or pre-fitted and then fixed because of sensitivity analysis results (cf. section 2.7). a negative feedback regulated by the HIF-1 complex (y 6 ). Hence, for the mTOR dynamics we obtained We denote phosphorylated STAT3 by y 8 . STAT3 basal activation and decay rates are a 8 and d 8 , respectively. Both IL-15 and mTOR are known to induce phosphorylation of STAT3 (29,34), here assumed to occur at rate k 6 and k 8 , respectively. All in all, the differential equation for STAT3 reads 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 7 ) and decay (d 7 ). NF-κB is further activated via IL-15 (27,28) (activation rate k 7 ), via mTOR (35) (k 15 ) and via the HIF-1 complex (36, 37) (k 14 ), yielding y ′ 7 (t) = a 7 + k 7 y 1 + k 14 y 6 + k 15 y 3 − d 7 y 7 .
HIF-1α mRNA basal synthesis and degradation are defined to occur at rate a 9 and d 9 , respectively. Further, we assumed that HIF-1α mRNA is regulated by NF-κB (at rate k 9 ) and STAT3 (at rate k 3 ), y ′ 9 (t) = a 9 + k 9 y 7 + k 3 y 8 − d 9 y 9 .
Following previous results (19), we assumed that asparaginyl hydroxylase FIH is at steady state (ϕ), whereas PHDs are In the right column we report the estimated parameter value, indicating mean and standard deviation (s.d.).
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 4 .
HIF-1β (y 5 ) is constitutively expressed by the cells (synthesis rate a 5 ), independently of the oxygen conditions (10). In hypoxia, HIF-1α accumulates and binds to HIF-1β (at rate k 4 ) forming the transcriptional complex HIF-1, which can dissociate (rate k 5 ). Hence, for HIF-1β and the HIF-1 complex we obtained the equations, and y ′ 6 (t) = k 4 y 4 y 5 − k 5 y 6 − d 6 y 6 , 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).  (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.

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 R 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 1 (0) = 1, whereas for unstimulated cells y 1 (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.

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, g(t, y(t), u, p).
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 ex j , y ex (t ex j ), u ex , p), i = 1, . . . , m, at measurement times t ex j ∈ [t 0 , t f ], j = 1, . . . , k ex , experimental data η ex ij are collected in each experiment ex. Experimental measurements contain additive noise where the errors ε ex ij are assumed to be statistically independent and normally distributed with zero mean value and variances (σ ex ij ) 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 0 ), hence to the initial condition y 0 . 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 0 and the model parameters must satisfy the steady state constraint,

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 2 -norm of the measurement errors, and further included a priori information p 0 by adding a Tikhonov regularization term (48) where the vector λ ∈ R 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

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 min s,p where the constraints F 2 (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 s k , p k+1 = p k + t k p k , t k ∈]0, 1] the increments s k , p k solve the linearized problem 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) 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.

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 0 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 R version 11.3 with the Ndsolve and manipulate routines, as well as DAESOL in VPLAN.

Sensitivity Analysis
In a second step, before starting parameter estimation, we performed a local sensitivity analysis at parameter values p 0 and corresponding solutions y ex (·|u ex , p 0 ), 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 0 = J(s 0 , p 0 ). If the sensitivity matrix J 0 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 0 . 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 0 . This allowed to fix 13 parameters which correspond to almost linear dependent columns of J 0 . The remaining 25 parameters are reliably identifiable and were estimated by mean of (PEP).

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 1 (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 2 , a 3 ) 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 R version 11.2.

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.
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.

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.

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 ≈ 10 −5 , k S ≈ 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 1 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, 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 1 model components in normoxia without (N, O 2 = 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 Experimentally collected data points are reported in Table S2.
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. 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.

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 1 in the interval [0, 10] nM h −1 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 1 = 0 nM h −1 ) and red curves to the highest value (a 1 = 10 nM h −1 ). 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 2 values in the presence of 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 3 (for mTOR) and a 8 (for STAT3) in the interval [0, 10] nM h −1 . 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 −1 , j = 3, 8) and red curves to the highest value (a j = 10 nM h −1 , 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.
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 3 ) and STAT3 (a 8 ) in the interval [0, 10] nM h −1 , 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 7 in the interval [0, 10] nM h −1 ) and STAT3 (a 8 in the interval [0, 10] nM h −1 ). 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 8 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 1 varying in [0, 10] nM h −1 ) and activation of STAT3 (a 8 varying in the interval [0, 10] nM h −1 ). The results confirm the synergistic effect of IL-15, STAT3 and DMOG in increasing HIF-1α levels.
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.  , 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 2 ± 25% or in the mTOR activation rate a 3 ± 25% 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).

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).
We compared the HIF-1α dynamics for the following     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).

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 nonlinear 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-15induced 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 2 ) 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).