Modeling the circadian regulation of the immune system: Sexually dimorphic effects of shift work

The circadian clock exerts significance influence on the immune system and disruption of circadian rhythms has been linked to inflammatory pathologies. Shift workers often experience circadian misalignment as their irregular work schedules disrupt the natural light-dark cycle, which in turn can cause serious health problems associated with alterations in genetic expressions of clock genes. In particular, shift work is associated with impairment in immune function, and those alterations are sex-specific. The goal of this study is to better understand the mechanisms that explain the weakened immune system in shift workers. To achieve that goal, we have constructed a mathematical model of the mammalian pulmonary circadian clock coupled to an acute inflammation model in the male and female rats. Shift work was simulated by an 8h-phase advance of the circadian system with sex-specific modulation of clock genes. The model reproduces the clock gene expression in the lung and the immune response to various doses of lipopolysaccharide (LPS). Under normal conditions, our model predicts that a host is more sensitive to LPS at circadian time (CT) CT12 versus CT0 due to a dynamic change of Interleukin 10 (IL-10), an anti-inflammatory cytokine. We identify REV-ERB as a key modulator of IL-10 activity throughout the circadian day. The model also predicts a reversal of the times of lowest and highest sensitivity to LPS, with males and females exhibiting an exaggerated response to LPS at CT0, which is countered by a blunted immune response at CT12. Overall, females produce fewer pro-inflammatory cytokines than males, but the extent of sequelae experienced by males and females varies across the circadian day. This model can serve as an essential component in an integrative model that will yield mechanistic understanding of how shift work-mediated circadian disruptions affect the inflammatory and other physiological responses.


The Sobol' sensitivity analysis
List of variables. We adopt the notation where names with a mix of upper and lower case letters (e.g., Per) denote mRNAs, and names in all caps (e.g., PER) denote proteins. Concentration of interleukin-6 IL10
Total number of activated phagocytic cells The initial condition is N (0) = 0.

Tissue damage marker
The initial condition is D(0) = 0.

Concentration of interleukin-10
The production of IL-10 in the basal state is represented by the constant sIL10. With N (0) = 0, the initial condition is IL- Tissue damage driven non-accessible interleukin-10 promoter The initial condition is Y IL10 (0) = 0.
Anti-inflammatory moderator At basal conditions, the system is assumed to be slightly anti-inflammatory. This was achieved by introducing a constant, sCA, into the ODE. Hence, with N (0) = 0, we obtain CA(0) = sCA dCA .

The Sobol' sensitivity analysis
The method of Sobol' [2] is a global and model independent sensitivity analysis method that is based on variance decomposition. Variance-based measures of sensitivity are attractive because they measure sensitivity across the whole input space, they can handle non-linear and non-monotonic functions and models. They can also measure the effect of interactions in such systems. A first order index, S i , is a measure for the variance contribution of the individual parameter to the total model variance, whereas a total order index, S T i is the result of the main effect of a given parameter and all its interactions with the other parameters. Note that for non-additive models as ours, interactions exist: S T i is greater than S i and the sum of all S i is less than 1. On the other hand, the sum of all S T i is greater than 1. By analyzing the difference between S T i and S i , one can determine the impact of the interactions between a parameter of interest and the other parameters [3]. In this work, we will not discuss how to implement the Sobol' method, we refer to [2,3] and the references therein.

The parameters considered for the sensitivity analysis
Based on our CJL male and CJL female models, 10 parameters are selected for the initial Sobol' sensitivity analysis. The set includes all parameters which took new values after fitting the CJL models. Lower and upper bounds to the parameter input space are obtained by halving and doubling nominal values, respectively (see Table M). Our outputs of interest are IL-6, TNF-α, IL-10 and D because they convey the most information about the state inflammation, and are also directly affected by the circadian clock model. For each output, we report first and total order effects. We performed 5000×(2+ number of input parameters) simulations to derive the Sobol' indices.

Results of the Sobol' sensitivity analysis
We calculated Sobol' indices in order to assess the relative influence of each parameter. These values measure the relative sensitivity of the outcome to each parameter (first-order) and to all the interactions with the other parameters (total-order). Fig A shows that total output uncertainty is primarily induced by the parameters kass pc, PER-CRY association rate, and vmax cry, Cry mRNA maximal transcription rate, and how they interact with the other parameters. The high sensitivity of D, IL-6, TNF-α and IL-10 to kass pc (Sobol index > 0.5), is due to the consequential role of the complex PER-CRY in the circadian circuitry. Inhibitor proteins PER and CRY dimerize to inhibit their own transcription as well as that of REV-ERB and ROR by acting on CLOCK-BMAL1 protein complex [4]. Overall, D, IL-6 and IL-10 are less sensitive with respect to the parameter set tested (see Table M), whereas TNF-α is more sensitive to parameters affecting CRY and its related complex, PER-CRY. Fig B shows the time course of the total-order Sobol' indices for vmax cry, kass pc and kass cb, the three most influtential parameters as shown in Fig A. D, IL-6 and IL-10 are sensitive kass pc and vmax cry throughout the inflammation, but the sensitivity of TNF-α to these parameters decreases over time. Moreover, the times of lower sensitivity of IL-10 to kass pc, correspond to the times of higher sensitivity of the cytokine to kass cb.
To assess the robustness of the outcomes to uncertainty in the input parameters in Table M, we plotted the computed simulation results showing the 90% quantile region for D, IL-6, TNF-α and IL-10 (see Fig C). Our results show that the inflammation output variables are qualitatively robust to uncertainty in the parameters that are modified by CJL.

Assessing sensitivities relative to coupling parameters
We extended our parameter input space to include all the clock model parameters and the coupling parameters. A sensitivity analysis of this larger parameter space allows us to assess how perturbations to the clock parameters can affect the output of the acute inflammation model, regardless of CJL. We performed 5000×(2+ number of input parameters) simulations to derive the Sobol' indices.
Our results indicate that D, IL-6, TNF-α and IL-10 are most sensitive to the coupling parameters, as should be expected. Fig D shows the Sobol' indices for the coupling parameters. These parameters directly link clock proteins to cytokines, and thus inferences about the effect of parameters on inflammation can be extended to inferences about the associated clock protein on inflammation.
TNF-α is more sensitive to xT N F CRY and xT N F ROR as CRY and ROR directly inhibit the production of the cytokine. IL-6 and IL-10 are sensitive to xIL6REV and xIL10REV , respectively. This is naturally explained by the direct inhibition of IL-6 and IL-10 by REV-ERB. We note that the damage marker, D, is sensitive to the same parameters as IL-6. This is to be expected because D is upregulated by IL-6. Fig E shows the time course of the total-order Sobol' indices for the coupling parameters. IL-6 and Damage. Both outcomes are sensitive to xIL6REV and xIL10REV throughout the duration of inflammation, but are also sensitive to xT N F CRY and xT N F ROR at the beginning of inflammation. This indicates that most of the effect of the clock on the damage marker is induced by REV-ERB, but there exists an initial joint effect of REV-ERB, CRY and ROR on IL-6 and D. CRY and ROR act indirectly through their modulation of TNF-α. TNF-α. The cytokine is particularly sensitive to xT N F CRY and xT N F ROR. The sensitivity of TNF-α to xT N F CRY decreases as its sensitivity to xT N F ROR increases during inflammation. Overall, ROR has a stronger influence on TNF-α compared to CRY. We note also the increased sensitivity of TNF-α to xIL10REV a few hours after the onset of inflammation. This is due to the inhibitory action of IL-10 on TNF-α. IL-10. Throughout the period of inflammation, IL-10 is most sensitive to xIL6REV and xIL10REV . The Sobol' index for xIL10REV decreases shortly after the onset of inflammation, before rising again five hours post-infection. This suggests a role for REV-ERB in the formation of the first peak in IL-10 as well as the second peak that occurs 5h after the beginning of inflammation (see Fig 5 in the manuscript). IL-10 is also sensitive to xT N F CRY and xT N F ROR at the start of inflammation. Since ROR and CRY inhibit TNF-α, the sensitivity of IL-10 to those parameters is because TNF-α peaks early during inflammation and activates its production.