Gulf War Illness: Is there lasting damage to the endocrine-immune circuitry?

Abstract

We reported previously that the persistence of complex immune, endocrine and neurological symptoms that afflict up to one third of veterans from the 1990-91 Gulf War might be supported by a misdirected regulatory drive. Here we use a detailed model of immune signaling in concert with an overarching circuit model of known sex and stress hormone co-regulation to explore how the failure of regulatory elements may further establish a self-perpetuating imbalance that closely resembles Gulf War Illness (GWI). Defects to the model were imparted iteratively and the stable regulatory modes supported by these altered immune-endocrine circuits were identified using repeated simulation experiments. In each case the predicted homeostatic regimes were compared to experimental data collected in male GWI (n=20 ) and matched healthy veterans (n=22 ). We found that alignment of GWI with a new homeostatic regime improved significantly when cortisol's normal anti-inflammatory activity was interrupted. Alignment improved further when this cortisol insensitivity was compounded by the loss of the normal antagonistic effects of Th1 cytokines on Th2 lymphocyte activation. Together these simulation results suggest altered glucocorticoid gene regulation compounded by possible changes in IGF-1 regulation of Th1:Th2 immune balance may be key underlying features of GWI.

Keywords:

• Gulf War Illness
• mathematical immunology
• numerical simulation
• regulatory control
• regulatory logic

Introduction

Gulf War Illness (GWI) is a complex multi-symptom illness associated with deployment to the Persian Gulf between 1990-91 and presenting with symptoms that manifest across several of the body's principal regulatory systems.^1-3 Our work and that of others strongly suggest that pathogenesis of this complex illness includes an immunologic ^4-7 and an endocrine component.⁸ Based on this, we proposed an illness model whereby GWI may involve a chronic imbalance in co-regulation between the nervous, endocrine and immune systems. To investigate this further we extended a discrete logic formalism first proposed by Mendoza and Xenarios⁹ to develop a coarse-grained but overarching model describing documented regulatory signaling within and between components of the hypothalamic-pituitary-adrenal (HPA), hypothalamic-pituitary-gonadal (HPG) and immune systems for both male and female physiology.¹⁰ This was followed by a more detailed fine grain model of the immune system incorporating oversight by stress hormone and male sex steroid signaling in a subsequent investigation.¹¹ In both cases logic circuits were assembled using an extension of representing normal regulatory physiology as documented in a broad survey of the medical literature.

Analysis of these logic circuit models revealed that the complexity associated with immune signaling and its oversight by sex and stress hormone regulation was such that these regulatory networks naturally supported more than one stable homeostatic regime. Comparison with a broad set of immune and endocrine markers measured in peripheral blood from GWI subjects^4,5 indicated that this chronic condition is proximal to such alternate homeostatic regimes and that homeostatic drive may be involved in the perpetuation of illness.^10,11 Specifically, we found that that the inclusion of basic immune function and male sex hormone regulation by the HPG axis with HPA function resulted in the emergence of 5 stable equilibrium states. Steroid and cytokine levels measured experimentally in healthy and ill veterans were compared to model predictions using a theoretical approximation of Fisher's method developed by Brown that estimates the overall significance of a set of p-values obtained from multiple independent tests of the same null hypothesis.¹² On this basis, male veterans with GWI showed poor alignment (p=0 .82) with the predicted nominal resting state, reinforcing that GWI represents a significant departure from normal endocrine-immune function. At this level of resolution, experimental data for GWI (p=0 .30) aligned best with a predicted state supporting chronically elevated cortisol levels, low testosterone and a shift toward Th1 immune activation. These results reinforced the premise that insults to the endocrine-immune system may displace normal homeostasis toward a more classical Th1 inflammatory signature with a concurrent and perhaps stabilizing endocrine component in GWI without requiring actual damage to the regulatory circuitry or the continued presence of the original insult. Nonetheless, while the immune and endocrine signatures observed in GWI aligned much more readily with these alternate stable regimes than with normal health this alignment was only partial. In other words while homeostatic drive may play a role in perpetuating GWI, the original insult may have also imparted some damage to the normal regulatory circuitry.

In this work we attempt to identify some of the affected processes by systematically simulating defects in regulatory circuits both at the broader endocrine-immune level¹⁰ as well as at the level of detailed immune signaling.¹¹ In each case the stable regulatory regimes supported by these altered circuits were identified and compared to available experimental data describing endocrine and immune marker expression in GWI. Significant improvement in alignment between alternate homeostatic states and GWI was obtained when the normal inhibition by cortisol of innate immune inflammatory activity was interrupted. Alignment improved further when this insensitivity to cortisol was compounded by a loss of responsiveness by Th2 lymphocytes to inhibition by Th1 cytokines resulting in a loss of Th1:Th2 polarization.

Methods

Literature Based Neuro-Endocrine Immune Signaling Networks. The Course Grain Model (CGM) used in this investigation is the male HPA-HPG Immune axis presented by Craddock et al.¹⁰ and has 14 biomarker nodes and 1 external influence (stress), with 40 interaction wires connecting them (Fig. 1). The Fine Grain Model (FGM), also for male HPA-HPG axis possesses a more detailed immune signaling network, has 21 biomarker nodes and 3 external influences (stress, infection, HPG drive), with 82 interaction wires connecting them (Fig. 2). Details of the FGM may also be found in previous work by our group.¹¹ The grouping of individual cytokines and hormones into the aggregate nodes is presented in Supplementary Tables S1, and S2 for the CGM and FGM respectively.

Figure 1. A coarse-grained model of endocrine-immune signaling. A coarse circuit logic representation of the principal signals linking the hypothalamic-pituitary-adrenal (HPA) and hypothalamic-pituitary-gonadal (HPG) endocrine axes with innate immune response (IIR) as well as T-helper 1(Th1) and 2 branches of the adaptive immune system as proposed by Craddock et al. (2014).¹⁰

Figure 2. A fine-grained model of immune cell signaling. A detailed circuit logic representation of immune cell signaling including the T-helper 1 (Th-1), 2, (Th2), 17, (Th17) as well as regulatory T-cell signaling (Treg) and their interface with sex and stress hormone regulation as proposed by Fritsch et al. (2013).¹¹

Discrete Logical Analysis. In this work we create a model of the immune system and its principal endocrine mediators by using a discrete ternary logical network analysis based on the analysis of Mendoza and Xenarios⁹ and Thomas,¹³ and developed in our previous analysis.^10,11 We encode the feedback mechanisms within the neuroendocrine-immune system based only on the direction (source and target) and type (activator or inhibitor) of interaction.

In this model, signaling molecules and cell types are represented as individual variables each capable of adopting 3 discrete state: −1 (down-regulated), 0 (nominal), and 1 (up-regulated). At any point in time t, the state of a system with N variables can be represented by the vector ${\displaystyle \overrightarrow{x}(t)}$, such that:(1) ${\displaystyle x(t)=(x_1(t),x_2(t),\dots ,x_N(t))}$(1) where ${\displaystyle x_i(t)}$ is the state of the i^th variable of the N variable system at time t. The image vector ${\displaystyle \overrightarrow{x}(t+1)}$ describes the preferred state toward which the system evolves in the next time increment. The state value of the image vector for the i^th variable is determined from its current state and a set of balanced ternary logic statements based on the current value of variable and the mode of action (i.e. activate or inhibit) of the neighboring input variables. These logic statements are expressed as follows (Eq. 2(2) ${\displaystyle x_i(t+1)=\{\begin{array}{l}(x_{i1}^A(t)\vee x_{i2}^A(t)\vee \mathrm{...}\vee x_{iX}^A(t))▿(x_{i1}^I(t)\vee x_{i2}^I(t)\vee \mathrm{...}\vee x_{iY}^I(t))\\ (x_{i1}^A(t)\vee x_{i2}^A(t)\vee \mathrm{...}\vee x_{iX}^A(t))\\ \neg (x_{i1}^I(t)\vee x_{i2}^I(t)\vee \mathrm{...}\vee x_{iY}^I(t))\end{array}}$(2) ):(2) ${\displaystyle x_i(t+1)=\{\begin{array}{l}(x_{i1}^A(t)\vee x_{i2}^A(t)\vee \mathrm{...}\vee x_{iX}^A(t))▿(x_{i1}^I(t)\vee x_{i2}^I(t)\vee \mathrm{...}\vee x_{iY}^I(t))\\ (x_{i1}^A(t)\vee x_{i2}^A(t)\vee \mathrm{...}\vee x_{iX}^A(t))\\ \neg (x_{i1}^I(t)\vee x_{i2}^I(t)\vee \mathrm{...}\vee x_{iY}^I(t))\end{array}}$(2) where the ∇, ∨, and ¬ symbols are ternary HIGH/LOW PASS, OR and NOT operators, ${\displaystyle x_{ij}^A}$ is the state of the i^th variable's j^th activator, ${\displaystyle x_{ik}^I}$is the state of the i^th variable's k^th inhibitor. The ternary operators given in Equation (2) are described in further detail in supplementary tables S3- S5. The first entry in Equation (2) is used when the variable possesses X activators and Y inhibitors, the middle when the variable has only X activators and last when the activator has only Y inhibitors.

Applying Equation (2) to each variable in the model for the m^th state of the system, ${\displaystyle \stackrel{\to m}{x}(t)}$, defines the image vector ${\displaystyle \stackrel{\to m}{x}(t)(t+1)}$ for the m^th state. With ${\displaystyle \stackrel{\to m}{x}(t)(t+1)}$ defined, the system may be updated asynchronously (allowing only one variable to change at a time) following the generalized logical analysis of Thomas.¹³ According to this method the i^th variable of the m^th state vector ${\displaystyle x_i^m(t)}$ is moved one step toward its preferred image ${\displaystyle x_i^m(t+1)}$ (e.g. If ${\displaystyle x_i^m(t)}$ = −1 and ${\displaystyle x_i^m(t+1)}$ = 1, then ${\displaystyle x_i^m(t+1)}$ is set to 0). Thus, for each current state of the system there are potentially several subsequent states toward which it may asynchronously evolve. Steady states are defined as those states for which the image vector is the same as the current state vector; in other words the state possesses an out degree of 0.

The number of molecular/cellular variables in the regulatory network model dictates the number of states the overall system is capable of possessing. Since each variable is capable of assuming one of 3 values (−1, 0, or 1), the number of states of the system is given by 3^N, where N is the number of biomarker variables. The CGM contains 15 soluble mediators and cell types, yielding 3¹⁴ = 4,782,969, (roughly 5 million) possible states, while the FGM has 3²¹ = 10,460,353,203 (roughly 10.5 billion) potential states. The above logical analysis is applied to each of these states to determine which states are stable.

GWI Cohort Sample Collection. Cytokine profiles and endocrine measures were obtained as part of a larger ongoing study in n = 27 GWI, 29 health control (HC) subjects recruited from the Miami Veterans Administration Medical Center. Subjects were all male with an average age of 43 y and BMI of 28. Inclusion criteria was derived from Fukuda et al.,¹⁴ and consisted in identifying veterans deployed to the theater of operations between August 8, 1990 and July 31, 1991, with one or more symptoms present after 6 months from at least 2 of the following: fatigue; mood and cognitive complaints; and musculoskeletal complaints. Subjects were in good health prior to 1990, and had no current exclusionary diagnoses.¹⁵ The use of the Fukuda definition in GWI is supported by Collins et al.¹⁶ Control subjects consisted of gulf war era sedentary veterans and were matched to illness subjects by age, body mass index (BMI) and ethnicity. Endocrine measures are comprised of plasma TEST and salivary CORT levels. Cytokine profiles consist of plasma levels of IL-1α, IL-1β, IL-2, IL-4, IL-5, IL-6, IL-8, IL-10, IL-12p70, IL-13, IL-17, IL-23, IFN-γ, TNF-α, and TNF-β. Summary statistics are presented in Supplementary Table S6. Additional details regarding this cohort and the laboratory assays performed are also available in Broderick et al.⁶

Statistical Comparison against Clinical Data Using Brown's Method. Brown's theoretical approximation¹² of Fisher's statistics was used to calculate the significance of alignment between experimental data and a given model predicted state. Fisher's method, a meta-analysis technique, combines probabilities to obtain the overall significance of a set of p-values obtained from independent tests of the same null hypothesis. The combined χ² statistic,(3) ${\displaystyle T_0=-2{\sum }_{i=1}^{N}ln(p_i)}$(3) where N is the number of measureable variables and p_i is the corresponding p-values under the null hypothesis, has a χ² distribution with 2N degrees of freedom assuming that the performed tests are independent. As the molecular variables of the endocrine and immune system interact with one another, as evidenced by the above connectivity diagrams, they are not independent. As a result, direct application of this test statistic is invalid, since the assumption of independence is violated. Brown¹² suggested a method for combining non-independent tests. If the tests are not independent, then the statistic T₀ has mean m = 2N and variance (σ²) given as,(4) ${\displaystyle {\sigma }^2=4N+2{\sum }_{i=1}^{N-1}{\sum }_{j=i+1}^{N}cov(-2ln{p}_i,-2ln{p}_j)}$(4) where p_i and p_j are the p-values for each test and the covariance (cov) is calculated as,(5) $cov(-2\text{ln}{p}_i,-2\text{ln}{p}_j)=\{\begin{array}{ll}{\rho }_{ij}(3.25+0.75{\rho }_{ij}),\hfill & 0\le {\rho }_{ij}\le 1\hfill \\ {\rho }_{ij}(3.27+0.71{\rho }_{ij}),\hfill & -0.5\le {\rho }_{ij}\le 0\hfill \end{array}$(5)

with ρ_ij being the unadulterated correlation between variable i and variable j. Finally, the overall significance P of a set of non-independent tests is calculated using the statistic T which under the null hypothesis follows the central χ² distribution, where T = T₀/c with 2N/c degrees of freedom and c = σ²/4N.

Here, we test if each experimental measure aligns with a given model predicted state. Our null hypothesis is that the experimental measures do not align. P-values for individual variables, p_i, are calculated using 2-sample t-tests between ill subjects and healthy controls. Where model predictions give a variable as high (+1), ‘right-handed’ one-tailed test are used, whereas a ‘left-handed’ test was used when model predictions are low (−1), to give the probability of obtaining the predicted value when the null hypothesis is true. For the case where the model predicts normal behavior for a variable (0) a 2-tailed t-test is used. However, the p-value from the 2-tailed test, p_two-tail, gives the probability that there is an observable difference between illness and control, which is the null hypothesis. To rectify this, when comparing to a model predicted variable of 0 we take the P-value to be p_i = 1 − p_two-tail, giving the probability of obtaining the predicted value when the null hypothesis is true. All cohort data was normalized using a Log2 transformation before T-tests and correlation calculations were performed. The unadulterated correlation values ρ_ij between 2 variables i and j were calculated in healthy subjects as the pairwise Pearson's linear correlation coefficient between variables. Where model variables represent an aggregate set of markers each experimentally measured constituent marker was compared individually to the model predicted value.

Systematic interruption of molecular signaling. Using Brown's method, as described above, a base proximity between GWI and the stables states predicted based on the complete CGM and FGM was established. Following this each complete model was used to create new models for proximal comparison to GWI first by deleting each individual connection in turn and then by deleting connections 2 at a time. These new altered networks were then analyzed using the discrete logical analysis, described above, to find alternate homeostatic states. Once all of the alternate states for the deletion models were found, they were compared to GWI via Brown's method,¹² and their subsequent p-value was noted.

A network is a set of n nodes, and m edges. To create new networks, the original set A was modified using the following algorithm:

For each edge i in A,

remove edge i,

save new set as set B,

return edge i.

Creating a new set B, of n nodes, and m-1 edges.

Then using the new set B,

For each edge i in B,

remove edge i,

save new set as set C,

End

Consider a complete set A with M edge elements G defined as,(6) ${\displaystyle A=A_M=\left\{G_1,G_2,\mathrm{...},G_M\right\}}$(6)

removal of a single edge interaction can then be defined as,(7) ${\displaystyle A_{M-i}=A_M-\left\{G_i\right\}=\left\{G_1,G_2,\mathrm{...},G_{i-1},G_{i+1},\mathrm{...},G_M\right\}}$(7) and a double deletion as(8) ${\displaystyle A_{M-i-j}=A_M-\left\{G_i,G_j\right\}=\left\{G_1,G_2,\mathrm{...},G_{i-1},G_{i+1},\mathrm{...},G_{j-1},G_{j+1},\mathrm{...},G_M\right\}}$(8)

and so forth.

Thus, we test the subsets of A_M defined as,(9) ${\displaystyle \begin{array}{lll}{A}_{M-i}& where& 1\le i\le M\end{array}}$(9) (10) ${\displaystyle \begin{array}{lll}{A}_{M-i-j}& where& 1\le (i,j)\le M\end{array}}$(10)

Ethics Statement

All subjects signed an informed consent approved by the Institutional Review Board of the University of Miami. Ethics review and approval for data analysis was also obtained by the IRB of the University of Alberta.

Results

Using Brown's combined null probability,¹² as described in the Methods section, a baseline statistical proximity between the observed GWI marker profile and the predicted stables states was calculated for a series of defects applied to (i) a broader but more coarse model of HPA-HPG-immune regulation (Fig. 1) as well as to (ii) a more detailed immune signaling network model (Fig. 2). Defects were applied individually and 2 at a time in order to systematically create a set of altered network models. In each case the stable states associated with this new altered wiring were identified using the discrete logical analysis described above and were compared for statistical similarity to GWI via Brown's method. Note that the reference state vectors corresponding to the healthy control and Gulf War illness conditions were based on the average immune and endocrine signatures shown in Supplemental Table S6.

Assessing potential defects in endocrine-immune co-regulation. Using the coarse but inclusive model of HPA-HPG-immune regulation¹⁰ (Fig. 1) a new series of endocrine-immune networks was created by deleting in turn each individual network connection. This resulted in the creation of 40 unique networks with a single interaction removed. This was followed by a second series of networks created by deleting one additional interaction from the newly generated networks. After removal of duplicate networks, this created 818 unique new networks, each with 2 connections removed. Recall that the original HPA-HPG-immune network model with all connections intact supported 4 alternate stable steady-states in the case of male physiology. As described in the Methods section stable states were identified exhaustively from all possible systems states supported by the network logic such that the system's next allowable state remains unchanged, or constitutes a state node with a zero out degree (Eq. 2(2) ${\displaystyle x_i(t+1)=\{\begin{array}{l}(x_{i1}^A(t)\vee x_{i2}^A(t)\vee \mathrm{...}\vee x_{iX}^A(t))▿(x_{i1}^I(t)\vee x_{i2}^I(t)\vee \mathrm{...}\vee x_{iY}^I(t))\\ (x_{i1}^A(t)\vee x_{i2}^A(t)\vee \mathrm{...}\vee x_{iX}^A(t))\\ \neg (x_{i1}^I(t)\vee x_{i2}^I(t)\vee \mathrm{...}\vee x_{iY}^I(t))\end{array}}$(2) ).

In the first alternate state (SS1) glucocorticoid receptor expression (GRD and GR) was elevated at low levels of ACTH, while in the second state (SS2) innate and Th1 immune responses were suppressed (low ICell, IIR, T1Cell, and T1Cyt) under elevated Th2 activity (high T2Cell and T2Cyt). Combining these features, the third state (SS3) was characterized by low ACTH at high glucocorticoid receptor expression (GR and GRD) and a shift toward Th2 immune activation (T2Cell and T2Cyt) with suppressed innate and Th1 activity (ICell, IIR, T1Cell and T1Cyt). The final state (SS4) best resembled GWI with a Brown's combined null probability of p=0 .30 as reported previously.¹⁰ This alternate homeostasis was characterized by hypercortisolism, low testosterone (TEST) levels and a shift toward a polarized Th1 immune response (low T2Cell, T2Cyt, GnRH, LH/FSH and TEST/EST, and high CORT, GRD, GR, T1Cyt and T1Cell). Of the 40 networks obtained by imparting a single deletion at a time to this model, one circuit supported an alternate stable state with a proximity to GWI corresponding to a combined Brown's p-value of 0.13 (Fig. 3) (see also SS1 in supplementary Table S7). This regulatory circuit was obtained by deleting the normal inhibitory actions of cortisol (Cort) on the inflammatory response supported by innate immune cell (Icell) activation. Proximity to GWI improved slightly when 2 interactions are removed yielding 2 altered circuits, each supporting a stable state with a combined Brown's p-value of divergence from GWI of p=0 .10 (Fig. 4). The first of these circuits is obtained by combining the above-mentioned defect in cortisol signaling with a loss of Th1 cytokine (T1Cyt) induced inhibition of Th2 cell (T2Cell) activity (Fig. 4(B) top panel). The second competing circuit corresponds to loss in cortisol anti-inflammatory signaling aggravated with a defect in Th2 cytokine (T2Cyt) production by the Th2 helper cell population (T2Cell) (Fig. 4(B) bottom panel).

Figure 3. Single defect in HPA-HPG-immune signaling. Statistical difference between observable marker levels in GWI and the alternate steady-state produced in the coarse grained HPA-HPG-Immune model with a single defect to the signaling network. Stable states were identified exhaustively from all possible systems states supported by the network logic such that the system's next allowable state remains unchanged (i.e., x(t+1)=x(t)), or constitutes a state node with a zero out degree (Eq. 2(2) ${\displaystyle x_i(t+1)=\{\begin{array}{l}(x_{i1}^A(t)\vee x_{i2}^A(t)\vee \mathrm{...}\vee x_{iX}^A(t))▿(x_{i1}^I(t)\vee x_{i2}^I(t)\vee \mathrm{...}\vee x_{iY}^I(t))\\ (x_{i1}^A(t)\vee x_{i2}^A(t)\vee \mathrm{...}\vee x_{iX}^A(t))\\ \neg (x_{i1}^I(t)\vee x_{i2}^I(t)\vee \mathrm{...}\vee x_{iY}^I(t))\end{array}}$(2) ). A single defect reduces Brown's combined null probability of a difference between GWI and the new modified model from a previous p-value of 0.30 to a value of 0.13 (panel A). This defect corresponds to interruption of normal cortisol-induced suppression of innate cell inflammatory activity (panel B).

Figure 4. Multiple defects in a coarse-grained signaling circuit. Interrupting two signals at a time in the coarse-grained HPA-HPG-immune network yields 2 solutions, both involving interruption of cortisol suppression of innate inflammatory response. In the first case this is accompanied by the loss of Th1 cytokine inhibition (T1Cyt) of Th2 cell activation (T2Cell) and in the second case by interruption of Th2 cytokine production (T2Cyt) by Th2 cells (T2Cell). As before, stable states were identified exhaustively from all possible systems states supported by the network logic such that the system's next allowable state remains unchanged, or constitutes a state node with a zero out degree. In both cases double interaction deletions produced stable states with a combined null p-value according to Brown's method of p=0 .10 for divergence with observed GWI marker levels.

Imparting defects in an intricate immune signaling network. The same procedure was performed on the fine-grained discrete logic model (FGM) of immune signaling presented previously in Fritsch et al.¹¹ (Fig. 2). Imparting individual defects created 82 new immune signaling networks. Interfering with 2 immune signals at a time produced 3327 unique networks. As before stable states were identified exhaustively from all possible systems states supported by the network logic such that the system's next allowable state remains unchanged (Eq. 2(2) ${\displaystyle x_i(t+1)=\{\begin{array}{l}(x_{i1}^A(t)\vee x_{i2}^A(t)\vee \mathrm{...}\vee x_{iX}^A(t))▿(x_{i1}^I(t)\vee x_{i2}^I(t)\vee \mathrm{...}\vee x_{iY}^I(t))\\ (x_{i1}^A(t)\vee x_{i2}^A(t)\vee \mathrm{...}\vee x_{iX}^A(t))\\ \neg (x_{i1}^I(t)\vee x_{i2}^I(t)\vee \mathrm{...}\vee x_{iY}^I(t))\end{array}}$(2) ). The original unaltered network supported 2 alternate stable attractors, one of which was separated from GWI by a combined Brown's null probability of p = 0 .13. The latter corresponded to a chronic state of Th17 immune activation in the context of elevated cortisol and low testosterone (TEST).

Of the 82 networks produced by interrupting individual immune signals 2 of these altered circuits supported alternate stable states that were more proximal to GWI (Fig. 5(A)). The first supported a stable state with a combined p-value for divergence from GWI of p=0 .10. The second circuit supported an alternate stable state with an even better alignment to GWI, namely a combined null probability of divergence of p=0 .05. The latter involved the loss of Th1 cytokine-induced (CK1) inhibition of Th2 helper cell activity (Th2) (Fig. 5(B)) (see also SS1 in supplementary Table S8). This solution reinforces the result obtained at coarser resolution with the HPA-HPG-immune multi-axis model. Removal of 2 connections at a time in this detailed immune model did not yield circuits supporting a better alignment with GWI (Figs. 6 and S1). Of the 3,321 dual-cut configurations created (82 x 81 x 1/2), 74 models supported 5 unique stable states that approximated GWI with a Brown's p=0 .05 (Fig. 6). All of these 74 circuit configurations included deletion of the inhibitory action of CK1 on the Th2 cells population. The remaining 7 of 81 available models containing a primary CK1/ Th2 deletion did not achieve the minimal Brown's p-value of 0.05. These 7 secondary defects are listed in Figure S1. Among these we find components of Th17 signaling. It is important to note that these findings do not indicate that Th17 signaling is up or downregulated but instead simply reinforce that these component mechanisms of normal Th17 signaling remain intact. Indeed their disruption decreases model fit to GWI.

Figure 5. Single compounding defect in detailed immune network circuitry. The statistical difference between observable marker levels in GWI and the alternate steady state produced with the introduction of a single defect in a detailed model of immune cell signaling. This defect reduces Brown's combined null probability of a difference between GWI and the new modified model from a previous p-value of 0.13 to a value of 0.05 (panel A). This defect corresponds to loss of Th1 cytokine (CK1) inhibition of Th2 cell activation (Th2) (panel B). Stable states were identified exhaustively from all possible systems states supported by the network logic such that the system's next allowable state remains unchanged, or constitutes a state node with a zero out degree.

Figure 6. Multiple defects in a detailed immune signaling circuit. Interrupting two signals at a time in the detailed immune network yields 5 stable state solutions. These do not improve the fit to the GWI chronic immune profile over that obtained with a single defect. Of the 3,321 possible double interaction deletions, 74 altered circuits supported stable states with a Brown's combined null p-value of p=0 .05 for divergence with observed GWI marker levels. As before stable states were identified exhaustively from all possible systems states supported by the network logic such that the system's next allowable state remains unchanged, or constitutes a state node with a zero out degree.

Network defects compromise robustness to challenge. In both of the damaged circuits mentioned above the defects modify the attractor landscape in a way that supports the emergence of attractors that more closely resemble GWI as a self-sustained chronic illness state. It is important to note that in each case even though a new attractor landscape emerges, this new landscape continues to support access to a stable healthy regulatory condition as shown in supplementary tables S7 and S8 (stable state SS3). This reinforces the notion that despite these alterations to the system some degree of stable remission may be possible. Indeed our simulation results suggest that it is not the outright accessibility of normal regulatory control that is compromised but rather the stability of this healthy state given these defects to the system. To explore this we simulated the response of the coarse grained model (Fig. 1) both with and without defects to idealized impulse challenges involving a sudden increase in stress, an acute immune response to injury and a combination of both stress and injury. Acute stress was simulated as a sudden impulse increase in cortisol (Cort) concentration whereas acute response to injury was simulated as a sudden elevation in innate immune signaling (IIR). Probabilistic discrete event simulations were repeated one million times in each case. Results presented in Supplementary Figure S2 show that both altered and unaltered circuits returned to healthy homeostasis roughly 2 out of 3 times in response to a sudden stressor. Once again both altered and unaltered circuits showed similar robustness to acute injury (impulse increase in IIR) with each returning to healthy homeostasis instead of escaping to a pathological steady state in roughly 1 out of 3 cases. However when simulated stress and injury were combined the unaltered HPA-HPG-immune circuit was much more robust to challenge, basically returning to a healthy control regime 1.5 times more frequently than in the case of the defective circuit. These results suggest that these defects may compromise the robustness of normal regulatory control and facilitating relapse but may not prohibit remission outright. It is important to note that the current models capture connectivity only and that a more thorough analysis of transition states separating attractors would require formal inclusion of response dynamics. Consequently these possible differences in local stability should be interpreted in relative terms only.

Discussion

In this work we continue to explore the involvement of homeostatic regulation in the perpetuation of symptoms in the chronic complex illness known as Gulf War Illness. In our earlier work we report discrete logic models of detailed immune signaling¹¹ as well as oversight of immune activation by stress and sex hormone regulation.¹⁰ In both cases logic circuits were assembled representing normal regulatory physiology as documented in a broad survey of the medical literature. As expected, we found that these regulatory networks naturally supported more than one stable homeostatic regime and that immune and endocrine signatures observed in GWI aligned much more readily with these altered stable regimes than with normal health. However this alignment was only partial suggesting that the original insult leading to GWI may have also imparted some lasting changes to the normal regulatory circuitry. To explore this we conducted an analysis of stable regulatory regimes supported when select defects were imparted to these circuit models. Using an exhaustive series of simulations based on a detailed model of immune signaling as well as a broader model of immune, stress and sex hormone regulation we found that interruption of cortisol's anti-inflammatory actions on innate immune activation produced an alternate regulatory program that more closely mimicked GWI. This similarity improved further when cortisol insensitivity was compounded with blunted inhibition of Th2 activity by the opposing Th1 immune axis.

Experimental data reported previously by this group⁶ indicate that GWI subjects often present with chronically elevated levels of cortisol (CORT). This same data also shows a loss of Th1:Th2 polarity as evidenced by concurrently elevated IL-13 (Th2) and TNF-β (Th1). Cortisol acts on naïve T-cells to repress both the Th1 promoting transcription factor, T-bet, and the Th2 promoting transcription factor, GATA-3. However, since the suppressive effect on T-bet is much stronger than the suppressive effect on GATA-3, a Th2 shift is ultimately favored¹⁷ under normal circumstances. The model presented here points to impaired glucocorticoid receptor (GR) signaling leading to glucocorticoid insensitivity and the persistence of elevated cortisol levels. Regulation of glucocorticoid responsiveness is complex however there is evidence that the latter is also reduced by exposure to pro-inflammatory cytokines, such as tumor necrosis factor α (TNF-α) and interleukin-1β (IL-1β),¹⁸ both of which are suggested to play a role in GWI.^4,19 Among its broad-reaching effects GR interacts with receptor tyrosine kinase to promote BDNF-induced release of glutamate. SiRNA transfection inducing a decrease in endogenous GR has been shown to decrease BDNF-dependent binding of phospholipase C-γ (PLC-γ) to tyrosine kinase B (TrkB) resulting in the suppression of neurotransmitter release. Consistent with the proposed GR insensitivity, our previous analysis of gene expression profiles in peripheral blood mononuclear cells (PBMC) indicated a decrease in Trk pathway activity in this cohort that were linked to changes in pathways supporting neurogenesis and axonal guidance.⁶ GR responsiveness is a known target of epigenetic control. Glucocorticoid resistance has been linked to reduced histone deacetylase-2 expression leading in turn to decreased repression of inflammatory genes.¹⁸ Moreover exposure to traumatic events has been shown to modify methylation of the glucocorticoid receptor gene NR3C1 promoter region with hypermethylation found in relation to early life stressors.²⁰ Conversely Yehuda et al.²¹ recently reported hypomethylation of NR3C1 promoter region in a cohort of male Gulf War veterans with PTSD, suggesting that epigenetic modification of the GR and associated changes in cortisol sensitivity the may depend on the timing, type and duration of exposure to adversity. Our simulation results would in principle support a hypermethylation of the glucocorticoid receptor gene promoter region leading to GR insensitivity and chronically elevated cortisol levels. This divergence in HPA status of GWI from PTSD is supported by Golier et al. (2007).⁸

Though our data shows a broad range of values for IFN-γ levels in peripheral blood from GWI subjects⁷ , this Th1 cytokine tends to be elevated in this illness group (p=0 .07 in a 2-way ANOVA).⁶ One would expect a trend toward lower expression of Th2 cytokines in GWI but these remain essentially normal in the case of IL-4 and IL-5, or increased in the case of IL-13, suggesting insensitivity to the normal IFN-γ induced down-regulation of Th2 polarization.²² In T helper cells down-modulation of IFN-γR2 acts as a negative regulatory mechanism that attenuates IFN-γ/signal transducer and activator of transcription-1 (STAT-1) signaling and limits the apoptotic effect of IFN-γ. This loss of IFN-γR2 that imparts resistance to IFN-γ is induced by IFN-γ itself in Th1 fated cells. This IFN-γR2 internalization is not observed in established Th2 clones. However other factors such as hormones may play a critical role in modulating IFN-γR2 expression and IFN-γ responsiveness in human T lymphocytes. For example IGF-1 is known to signal for IFN-γR2 internalization and limit IFN-γ/STAT-1 signaling in human T cells.²³ Though we did not measure IGF-1 levels in our previous work, we did find high IL-13 (p<0.01) which is known to increase IGF-1 production in macrophages.²⁴ We also found genomic evidence of increased NF-kB activation in this GWI cohort which promoted by IGF-1 via phosphatidylinositol 3-kinase (PI3K) / Protein Kinase C (PkC) engagement.²⁵ While not conclusive these simulation results suggest a mechanism by which an IGF-1 rich environment may abrogate IFN-γ / STAT1 mediated apoptosis of Th2 cells leading to a joint Th1/Th2 profile observed experimentally. Moreover IGF-1 has been found to inhibit whereas IGF-2 potentiated K+-evoked ACh release from hippocampal slices in rat.^26,27 A leading hypothesis in GWI involves extended low-level exposure to Acetylcholine esterase inhibitors²⁸ and elevated IGF-1 levels might be expected response to high ACh levels. This modulating action of IGF-1 was not incorporated a priori in the model and may explain why the alignment to GWI is enhanced if we provide a means for sustained Th2 activity in an environment trending toward increased IFN-γ and TNF-β.

These simulation results suggest that transient insult to the endocrine-immune regulatory network in GWI, in addition to prompting a persistent adaptive response, may have left lasting effects, perhaps through epigenomic mechanisms, to GR expression that impair HPA axis and Th1/Th2 balance. The latter may also be mediated by chronic increases in IL-13 induced IGF-1 expression. The latter response may have been prompted initially as a response to elevated ACh if exposure to AChE inhibitors was the insulting agent. In addition to suggesting possible pathogenic processes in GWI, these simulation results inform future improvements in the computational model in particular direct and preferential suppression of Th1 by GC should be modeled in greater detail, as should the immuno-modulatory effects of IGF and their relation to ACh expression. These efforts are ongoing.

It is also interesting to note that while these defects produced stable attractors presenting with immune and endocrine profiles that more closely resembled GWI, the resulting landscape still supported access to a normal homeostatic regime. Exploratory simulations suggested however that this healthy regulatory state was somewhat less robust to insult and that the defects imparted to the networks made it easier to escape health and lapse into a persistent pathological state. This result is not all that surprising when we consider the inherent properties of biological networks, which select simultaneously for adaptability and robustness.²⁹ Defects such as these might be compared to non-lethal deletions in gene regulatory networks. While the current connectivity models strictly enforce the sequence of allowable states they do not account for transition dynamics. While they remain exploratory these early findings nonetheless motivate the inclusion of more detailed transition dynamics in a future iteration of these models in order to support of a more thorough analysis of the topology in regions that lie between stable attractor states.

Disclosure of Potential Conflicts of Interest

The authors declare that they have no conflicts of interest.

Supplemental Material

Supplemental data for this article can be accessed on the publisher's website.

Author Contributions

MAR, RMdR, TJAC, JZ and GB performed computational simulations and analysis of computational results. TJAC, NGK, MAF and GB, compiled physiological data and developed the immune-endocrine connectivity network. MAR, TJAC and GB performed comparisons between experimental and computational data and interpreted the results. MAF, NGK and ZB collected and analyzed the experimental data. TJAC, MAR, RMdR, and GB conceived of the discrete logical analysis algorithm. All authors contributed to the drafting of this article, and approve of its contents.

Acknowledgments

This research was conducted in collaboration with the High Performance Computing team at the University of Miami Center for Computational Science (CCS) (http://ccs.miami.edu). The authors would also like to thank Dr. Andrée Gruslin for many helpful discussions regarding the many regulatory effects of insulin growth factor (IGF).

Funding

Funding was provided by US Department of Defense Congressionally Directed Medical Research Program (CDMRP) awards (http://cdmrp.army.mil/) GW093042 (Broderick - PI) and GW080152 (Klimas - PI) as well as Merit Awards from the US Veterans Affairs (Klimas - PI).