Patent Publication Number: US-2022230703-A1

Title: Predicting clinical intervention strategy for treatments of a complex disease

Description:
CROSS REFERENCE TO RELATED APPLICATIONS 
     This application claims priority to PCT Patent Application No. PCT/US2020/035367, filed on May 29, 2020, entitled “Computational System And Process For Predicting A Clinical Intervention Strategy For Treatment Of A Complex Disease,” which claims priority to U.S. Provisional Patent Application No. 62/854,196, filed on May 29, 2019, entitled “Computational Pipeline And Method For Prediction of Clinical Intervention Strategies For Complex Diseases”, the teaching of each are hereby incorporated by reference in their entirety. 
    
    
     FEDERALLY SPONSORED RESEARCH AND DEVELOPMENT 
     This invention was made with government support under W81XWH-10-1-0774, W81XWH-13-2-085, W81XWH-16-1-0552, W81XWH-15-1-0582, and W81XWH-16-1-0632 awarded by the Department of Defense and NS090200 awarded by the National Institutes of Health. The government has certain rights in the invention. 
    
    
     FIELD OF THE INVENTION 
     The invention is encompassed within the fields of computational biology, network biology, and medicine and particularly relates to the use of advanced computing capabilities to design models and simulations of physical phenomena that aid in understanding and solving complex biological problems and most particularly relates to pipeline computing applications that integrate multiple types of scientific analyses to enable prediction of clinical intervention strategies useful in treatment of complex and/or multifactorial diseases. 
     BACKGROUND 
     The present invention generally relates to machine learning and more specifically relates to predictive modeling of treatment courses for a complex disease designed to return the overall neuroimmune system to healthy regulation. 
     Gulf War Illness (GWI) is a chronic disorder that affects up to one-third of the 700,000 veterans of the Persian Gulf War of 1991. Many symptoms of GWI, including fatigue, musculoskeletal pain, cognitive, and neurological problems [ 1 ], are indicative of a neuroinflammatory disorder, and abnormal regulation of the hypothalamic-pituitary-adrenal (HPA) axis, a key component in the body&#39;s stress response, has been associated with this chronic condition [2-6]. Yet, the cause and underlying mechanisms of GWI are still relatively unknown, and there is no known treatment. Progress needs to be made to provide relief to this aging population, but the heterogeneity of the illness presentation provides complications. Subtyping the illness is one option to improve diagnostics and identify potential subtype-specific treatment avenues. 
     Possible explanations for the cause of GWI include exposure to mild traumatic brain injury and chemical/biological weapons [7] as well as the involvement of a neuroinflammatory signaling cascade triggered by exposure to a neurotoxin (e.g., the organophosphatesarin) and exposure to a combat environment [8-10]. Genetic susceptibility [11,12] of veterans may further contribute to the resulting changes in brain function observed in GWI [13-15]. Furthermore, confounding and/or exacerbating factors often complicate this diagnosis. Veterans suffering from GWI are often also diagnosed with post-traumatic stress disorder (PTSD) [16,17]. A population-based study estimates that 35% of veterans with GWI also suffer from PTSD [18], which exacerbates all the overall clinical symptoms of GWI [19]. As there is a substantial overlap of PTSD symptoms with the symptoms defining GWI [20], with both PTSD [21-24] and GWI [6,25,26] exhibiting immunological dysfunction and abnormal neural synchrony [14], this suggests that GWI and PTSD share several core pathophysiological processes and that these processes may be altered in the context of co-morbidity. As such, the presence or absence of PTSD in GWI provides one mean to subtype this illness. 
     While the immune system, as an independent agent, provides the human body with a defense system that protects against foreign proteins and chemicals, it cannot be considered independent from the body as a whole or the brain in particular. The brain plays a critical role in the regulation and modulation of immunity, effectively integrating the nervous and immune systems, the two most complex systems of the body for the maintenance of homeostasis, into one. The study of these neuroimmune interactions along with the endocrine system and resulting behaviors is collectively known as psychoneuroimmunology and studying the functional significance of the relationships among these systems is the focus of the field. Analyzing the mechanisms that functionally alter this single overarching system can lead to a significantly enhanced understanding of the unique and combined contributions of hormones, the brain, behavior, and the immune system to health and disease. 
     The immune system and the brain communicate through several signaling pathways. Two of the major pathways involved in this crosstalk are the sympathetic nervous system via the sympathetic-adrenal-medullary axis and the hypothalamic-pituitary-adrenal (HPA) axis. Both release bloodborne mediators, including cytokines, chemokines, and hormones, which are then subject to further regulation. The brain and spinal cord, which comprise the central nervous system (CNS), are kept isolated from signaling in the blood flow of the peripheral body by a monocellular interface called the blood-brain barrier (BBB). This division is made possible by several parallel barriers that include most notably the capillary bed of the CNS and the choroid plexus. These barriers at one level create the dichotomy between the circulating factors of the immune system and the components of the CNS only to regulate the interactions between the immune system and the CNS at other levels. The BBB is thus an integral part of the neuroimmune axis. Many factors, including stress, regulate the permeability of the BBB. Corticotropin-releasing hormone (CRH), released by the HPA axis during times of stress, has proinflammatory effects mediated through mast cells, which regulate BBB permeability [27,28]. Furthermore, CRH release leads to the secretion of catecholamines and glucocorticoids, which exert their effects on the BBB. Exposure to acute or chronic stress activates the HPA axis, whether the stressor is mental or physiological distress, even inducing, and can lead to BBB breakdown [27,28]. 
     Not only is the HPA axis activated by stress, HPA activity is also tied to both the hypothalamic-pituitary-gonadal (HPG) axis and the immune system, among others, and is further influenced by brain signaling outside the hypothalamus. Moreover, weakening the BBB can lead to additional or enhanced interactions among these and other systems. This level of integrated connectivity is core to the body&#39;s adaptive capability and contributes to its resilience and vulnerabilities. As the HPA axis plays a crucial role in the body&#39;s response to stress [29-31] and is one of the major information pathways connecting the peripheral systems to the CNS, several mathematical models of HPA dynamics have been proposed to understand its function [32-37]. Few models account for its detailed feed-forward regulation [37], but this is highly relevant as signaling networks such as the HPA axis do not work in isolation. Even fewer account for its connection to other bodily systems. 
     A complex network of regulatory interactions that possess feedforward and feedback connections, such as the aforementioned neuroimmune axis, can have multiple response dynamics that may include several stable homeostatic states beyond a normal healthy state [38-46]. Due to homeostatic feedback mechanisms, acute or minor changes to this system return to the baseline, or a normal healthy state. Nevertheless, under more severe perturbations, the system may shift from a normal homeostatic state to a new stable state. 
     This disturbance may consist of a single event or an ongoing external insult that further changes the underlying dynamics. We believe that these regulatory alterations are a key cause of the resistance to therapy seen in neuroinflammatory illnesses and may explain some of the hardships in attaining lasting remission. By identifying and detailing these alternative homeostatic regimes, it may be possible to design treatments aimed at returning the overall system back to its original healthy state. 
     SUMMARY 
     The present invention builds on the inventors&#39; ongoing research directed at mapping complex inflammatory mechanisms in GWI to improve our understanding of the immunologic underpinnings of GWI, the compounding effects of co-morbidity with PTSD, and the potential of this co-morbidity to define a unique subtype of GWI. Furthermore, we propose to take advantage of our ongoing work in predictive modeling to assess possible changes to putative treatments of GWI in the context of a probable PTSD diagnosis or lack thereof. We do this by constructing a logic model representation of the neuroimmune system across the BBB, connecting the CNS to blood-based system. Examples of blood-based systems include is at least one of the subject&#39;s immune system, the subject&#39;s endocrine system, the subject&#39;s hormone system, or a combination thereof. As we cannot assess the neuroimmune state in the brain of living GWI subjects, we topologically compare the model results of a GWI cohort with and without PTSD symptoms to mouse models of GWI in blood and use the corresponding neuroimmune states in a mouse brain to infer the blood-brain state of the GWI subjects. This inference is then used to guide predictive modeling of treatment courses designed to return the overall neuroimmune system to healthy regulation. The present invention is also applicable to other complex diseases such as Parkinson&#39;s disease. 
     Gulf War Illness (GWI) is a persistent chronic neuroinflammatory illness exacerbated by external stressors and characterized by fatigue, musculoskeletal pain, cognitive, and neurological problems linked to underlying immunological dysfunction for which there is no known treatment. As the immune system and the brain communicate through several signaling pathways, including the hypothalamic-pituitary-adrenal (HPA) axis, it underlies many of the behavioral and physiological responses to stressors via bloodborne mediators, such as cytokines, chemokines, and hormones. Signaling by these molecules is mediated by the semipermeable blood-brain barrier (BBB) made up of a monocellular layer forming an integral part of the neuroimmune axis. BBB permeability can be altered and even diminished by both external factors (e.g., chemical agents) and internal conditions (e.g., acute or chronic stress, or cross-signaling from the hypothalamic-pituitary-gonadal (HPG) axis). Such a complex network of regulatory interactions that possess feedforward and feedback connections can have multiple response dynamics that may include several stable homeostatic states beyond normal health. Here we compare immune and hormone measures in the blood of human clinical samples and mouse models of Gulf War Illness (GWI) subtyped by exposure to traumatic stress for subtyping this complex illness. We do this via constructing a detailed logic model of HPA-HPG-Immune regulatory behavior that also considers signaling pathways across the BBB to neuronal-glial interactions within the brain. We apply conditional interactions to model the effects of changes in BBB permeability. Several stable states are identified in the system beyond typical health. Following alignment of the human and mouse blood profiles in the context of the model, mouse brain sample measures were used to infer the neuroinflammatory state in human GWI and perform treatment simulations using a genetic algorithm to optimize the Monte Carlo simulations of the putative treatment strategies aimed at returning the ill system back to health. We identify several ideal multi-intervention strategies and potential drug candidates that may be used to treat chronic neuroinflammation in GWI. 
     More specifically, the present invention provides a system, a computer program product, and a method for predicting a clinical intervention strategy to treat a complex disease such as Gulf War Disease or Parkinson&#39;s disease. The method begins accessing a set of predefined conditions that include a number of inhibitors, a number of activators, a healthy control value, and a conditional state of variables. Next, a signaling network is accessed that represents a subject across a blood-brain barrier connecting to the subject&#39;s central nervous system to a blood-based system of the subject. In one example, the blood-based system of the subject is the subject&#39;s immune system, the subject&#39;s endocrine system, the subject&#39;s hormone system, or a combination thereof. 
     The signaling network includes the conditional state of one or more variables that produce outputs based on inputs and activations due to external influences on the signaling network. Also, the signaling network includes a set of nodes interconnected by at least one of a static stimulatory edge, a conditional stimulatory edge, an inhibitory edge, or a combination thereof. In one example, the signaling network is a neuroimmune signaling network the further includes using the nodes in the set of nodes to model one or more of immune signaling molecules and immune cells, hypothalamic-pituitary-adrenal (HPA) axis and hypothalamic-pituitary-gonadal (HPG) axis, andneuronal-glia interactions, growth factors, blood-brain barrier, neurotransmitter, and neuroimmune signaling molecules. 
     A ternary logic modeling is accessed in which signaling molecules and cell types are represented as individual variables, each capable of adopting one of three distinct states. Next, a computational algorithm is performed that relies on repeated random sampling to obtain numerical results representing of state evolution of the signaling network to identify distinct states as potential stable behaviors compared with a healthy control. Next, the model results of the neuroimmune signaling network are topologically compared with and without behavioral disorder to mouse models of a complex disease in blood and use corresponding states in a mouse brain to infer the blood-brain state of the subjects. Based based on the conditional state of the variables equal to a given state, updating at least one of the conditional stimulatory edge to a static stimulatory edge and repeating step b through step e until convergence or for a set number of iterations. The inferred blood-brain state of the subject to provide treatment courses thereto. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Office upon request and payment of the necessary fee. 
       In the accompanying figures, like reference numerals refer to identical or functionally similar elements throughout the separate views. The accompanying figures, together with the detailed description below, are incorporated in and form part of the specification and serve to further illustrate various embodiments and to explain various principles and advantages all in accordance with the present invention, in which: 
         FIG. 1  is a Table 1 of the statistical comparisons were made between both the total GWI and HC groups (p 2 ) as well as between the GWI H , GWI L , and HC groups (p 3 ), using ANOVA for the continuous variables and the χ 2  test for categorical variables, according to one aspect of the present invention. 
         FIG. 2  is a graph illustrating human blood cytokine levels in PTSD-based GWI subtypes as compared to HC at rest (T0) and peak exercise during a maximal graded exercise challenge (T1), according to one aspect of the present invention. For each cytokine, the data is presented, from left to right as follows: GWIL T0; GWIH T0; GWIL T1; and GWIH T1. 
         FIG. 3A  through  FIG. 3C  are graphs showing cytokine expression in a long-term mouse model of Gulf War Illness, according to one aspect of the present invention. For  FIGS. 3B and 3C , the data is presented from left to right as set forth in the respective legends from top to bottom. 
         FIG. 4  is a theoretical male neuroimmune signaling network, according to one aspect of the present invention. 
         FIG. 5  illustrates predicted stable states of the theoretical male neuroimmune signaling network presented in  FIG. 4 , according to one aspect of the present invention. 
         FIG. 6A  through  FIG. 6D  are graphs showing overall topology of the stable state locations is retained between humans and mice, according to one aspect of the present invention. 
         FIG. 7A  through  FIG. 7D  are graphs showing overall topology of the stable state locations is retained between humans and mice, according to one aspect of the present invention. 
         FIG. 8A  through  FIG. 8D  is a comparison of the GWI trauma subtype blood profile topologies with mouse models under challenge, according to one aspect of the present invention. 
         FIG. 9A  through  FIG. 9E  are a comparison of the mouse model conditions brain and blood profile topologies, according to one aspect of the present invention. 
         FIG. 10  is a Table 2 illustrating a SS5 start state with a minimum of three treatment targets needed to obtain at least a 1% return to health, according to one aspect of the present invention. 
         FIG. 11  is a Table 3 illustrating a SS4 start state, again, a minimum of three treatment targets were needed to obtain at least 1% return to health, according to one aspect of the present invention. 
         FIG. 12  is a Table 4 illustrating a SS2 start state, only two interventions were needed to obtain at least a 1% return to health, according to one aspect of the present invention. 
         FIG. 13  illustrates a block diagram of a system for implementing the processes described above, according to one aspect of the present invention. 
     
    
    
     DETAILED DESCRIPTION 
     This application includes references denoted in brackets with year of publication. These references are listed in the information disclosure statement (IDS) filed herewith. The teachings of each of these references is hereby incorporated hereinto in their entirety. 
     Non-Limiting Definitions 
     The term “artificial neural network” or “ANN” is a learning system modeled after the human brain, operated by a large number of processors operating in parallel. 
     The term “hyperparameters” refers to parameters that define properties of the training model, but cannot be learned directly from the process of training the model. Hyperparameters are usually set before the actual training process begins and describe properties such as: the depth of a tree, the rate of learning, the number of hidden layers, or the number of clusters. They are also known as “meta parameters.” 
     The term “model parameters” refers to the parameters in a machine learning model. Model parameters are learned from training data. 
     The term “meta parameters” is another term for “hyperparameters.” 
     The term “training parameters” is another term for model parameters. 
     Introduction 
     Gulf War Illness (GWI) is a chronic disorder that affects up to one-third of the 700,000 veterans of the Persian Gulf War of 1991. Many symptoms of GWI, including fatigue, musculoskeletal pain, cognitive, and neurological problems [1], are indicative of a neuroinflammatory disorder and abnormal regulation of the hypothalamic-pituitary-adrenal (HPA) axis, a key component in the body&#39;s stress response, has been associated with this chronic condition [2-6]. Yet, the cause and underlying mechanisms of GWI are still relatively unknown, and there is no known treatment. Progress needs to be made to provide relief to this aging population, but the heterogeneity of the illness presentation provides complications. Subtyping the illness is one option to improve diagnostics and identify potential subtype-specific treatment avenues. 
     Possible explanations for the cause of GWI include exposure to mild traumatic brain injury and chemical/biological weapons [7], as well as the involvement of a neuroinflammatory signaling cascade triggered by exposure to a neurotoxin (e.g., the organophosphate sarin) and exposure to a combat environment [8-10]. Genetic susceptibility [11,12] of veterans may further contribute to the resulting changes in brain function observed in GWI [13-15]. Furthermore, confounding and/or exacerbating factors often complicate this diagnosis. Veterans suffering from GWI are often also diagnosed with post-traumatic stress disorder (PTSD) [16,17]. A population-based study estimates that 35% of veterans with GWI also suffer from PTSD [18], which exacerbates all the overall clinical symptoms of GWI [19]. As there is a substantial overlap of PTSD symptoms with the symptoms defining GWI [20], with both PTSD [21-24] and GWI [6,25,26] exhibiting immunological dysfunction and abnormal neural synchrony [14], this suggests that GWI and PTSD share several core pathophysiological processes and that these processes may be altered in the context of co-morbidity. As such, the presence or absence of PTSD in GWI provides one means to subtype this illness. 
     While the immune system, as an independent agent, provides the human body with a defense system that protects against foreign proteins and chemicals, it cannot be considered independent from the body as a whole, or the brain in particular. The brain plays a critical role in the regulation and modulation of immunity, effectively integrating the nervous and immune systems, the two most complex systems of the body for the maintenance of homeostasis, into one. The study of these neuroimmune interactions along with the endocrine system and resulting behaviors is collectively known as psychoneuroimmunology and studying the functional significance of the relationships among these systems is the focus of the field. Analyzing the mechanisms that functionally alter this single overarching system can lead to a significantly enhanced understanding of the unique and combined contributions of hormones, the brain, behavior, and the immune system to health and disease. 
     The immune system and the brain communicate through several signaling path-ways. Two of the major pathways involved in this crosstalk are the sympathetic nervous system via the sympathetic-adrenal-medullary axis and the hypothalamic-pituitary-adrenal (HPA) axis. Both release bloodborne mediators, including cytokines, chemokines, and hormones, which are then subject to further regulation. The brain and spinal cord, which comprise the central nervous system (CNS), are kept isolated from signaling in the blood flow of the peripheral body by a monocellular interface called the blood-brain barrier (BBB). This division is made possible by several parallel barriers that include most notably the capillary bed of the CNS and the choroid plexus. These barriers at one level create the dichotomy between the circulating factors of the immune system and the components of the CNS only to regulate the interactions between the immune system and the CNS at other levels. The BBB is thus an integral part of the neuroimmune axis. Many fac-tors, including stress, regulate the permeability of the BBB. Corticotropin-releasing hormone (CRH), released by the HPA axis during times of stress, has proinflammatory effects mediated through mast cells, which regulate BBB permeability [27,28]. Furthermore, CRH release leads to the secretion of catecholamines and glucocorticoids, which exert their effects on the BBB. Exposure to acute or chronic stress activates the HPA axis, whether the stressor is mental or physiological distress, even inducing, and can lead to BBB break-down [27,28]. 
     Not only is the HPA axis activated by stress, HPA activity is also tied to both the hypothalamic-pituitary-gonadal (HPG) axis and the immune system, among others, and is further influenced by brain signaling outside the hypothalamus. Moreover, weakening the BBB can lead to additional or enhanced interactions among these and other systems. This level of integrated connectivity is core to the body&#39;s adaptive capability and contributes to its resilience and vulnerabilities. As the HPA axis plays a crucial role in the body&#39;s response to stress [29-31] and is one of the major information pathways connecting the peripheral systems to the CNS, several mathematical models of HPA dynamics have been proposed to understand its function [32-37]. Few models account for its detailed feed-forward regulation [37], but this is highly relevant as signaling networks such as the HPA axis do not work in isolation. Even fewer account for its connection to other bodily systems. 
     A complex network of regulatory interactions that possess feed-forward and feed-back connections, such as the aforementioned neuro-immune axis, can have multiple response dynamics that may include several stable homeostatic states beyond a normal healthy state [38-46]. Due to homeostatic feedback mechanisms, acute or minor changes to this system return to the baseline or a normal healthy state. Nevertheless, under more severe perturbations, the system may shift from a normal homeostatic state to a new stable state. This disturbance may consist of a single event or an ongoing external insult that further changes the underlying dynamics. We believe that these regulatory alterations are a key cause of the resistance to therapy seen in neuroinflammatory illnesses and may ex-plain some of the hardships in attaining lasting remission. By identifying and detailing these alternative homeostatic regimes, it may be possible to design treatments aimed at returning the overall system back to its original healthy state. 
     In this work, we propose to build on our ongoing research directed at mapping com-plex inflammatory mechanisms in GWI to improve our understanding of the immuno-logic underpinnings of GWI, the compounding effects of co-morbidity with PTSD, and the potential of this co-morbidity to define a unique subtype of GWI. Furthermore, we propose to take advantage of our ongoing work in predictive modeling to assess possible changes to putative treatments of GWI in the context of a probable PTSD diagnosis or lack thereof. We do this by constructing a logic model representation of the neuroimmune system across the BBB connecting the CNS to to blood-based system. Examples of blood-based systems include is at least one of the subject&#39;s immune system, the subject&#39;s endocrine system, the subject&#39;s hormone system, or a combination thereof. As we cannot assess the neuroimmune state in the brain of living GWI subjects, we topologically compare the model results of a GWI cohort with and without PTSD symptoms to mouse models of GWI in blood and use the corresponding neuroimmune states in a mouse brain to infer the blood-brain state of the GWI subjects. This inference is then used to guide predictive modeling of treatment courses designed to return the overall neuroimmune system to healthy regulation. 
     Results 
     1. The GWI Cohort 
     The total cohort used in this study was comprised of 89 male veterans recruited through the Miami Veterans Affairs (VA) Medical Health Center. The average age of the cohort was 43.68±7.42 years, with an average BMI of 29.49±5.42. The sample consisted of 2.2% Asian, 22.5% Black, 46.1% White Hispanic, 28.1% White, and 1.1% Other. The full demographics are presented in Table 1. 
     The full sample included n=51 male veterans meeting the criteria for GWI in addition to n=38 healthy sedentary veteran controls (HC). McDonald et al. (2014) [47] examined the diagnostic accuracy of applying a cut score method to the Davidson Trauma Score (DTS) in a population of U.S. military veterans and found that using a cut score of 70 for the total DTS score offered optimal diagnostic accuracy, correctly classifying 90% of cases and providing an accurate estimate of PTSD population prevalence. Of the 51 GWI subjects, n=32 had a total DTS of 70 or greater. The group of GWI subjects with DTS scores above 70 are considered to have a high probability of having PTSD and are denoted as GWIH, while the remainder have a low probability of having PTSD and are labelled as GWIL. 
     Turning to  FIG. 1 , shown in Table 1 are the statistical comparisons were made between both the total GWI and HC groups (p 2 ) as well as between the GWI H , GWI L , and HC groups (p 3 ), using ANOVA for the continuous variables and the χ 2  test for categorical variables. No statistical differences were found in age, BMI, racial representation, marital status, employment status, or an average number of years in school. 
     Blood plasma cytokine levels for the GWI subgroups at rest (T0) and peak effort in a maximal graded exercise test (T1) as compared to HC at rest are presented in  FIG. 2 . Statistical comparisons were made using heteroscedastic two-tailed t-tests between the GWI subgroups and HC at both timepoints, and between the GWI subgroups within each timepoint. While there appears to be trends of immune dysregulation in the GWI subgroups at rest, there were no statistical differences between the GWI subgroups and HC or between the two subgroup conditions, owing to the large degree of variability within the measures. After immune stimulation via the graded exercise challenge, the immune profiles of the GWI subgroups again show trends of immune dysregulation with several statistically significant differences, namely, elevated IL-8, IFN-γ, and TNFα in both subgroups, depressed IL-4 in GWI L , and elevated IL-1β in GWI H . Again, owing to the large degree of variability in the measures, no significant differences were found between the GWI subgroups at T1. 
     Turning to  FIG. 2  is illustrates human blood cytokine levels in PTSD-based GWI subtypes as compared to HC at rest (T0) and peak exercise during a maximal graded exercise challenge (T1). Individual cytokine values are normalized to the HC values, and bars represent the mean ratios±SEM for each cytokine. Significance was set at * p&lt;0.05, ** p&lt;0.01, and *** p&lt;0.001 as compared to HC. The false discovery rate, as calculated via Storey et al. [48], was &lt;0.1 for all significant measures. Values are found in Table 1 of  FIG. 1 . 
     2. A Long-Term Mouse Model of Gulf War Illness 
     Veterans with GWI had the potential to be exposed to several toxicants during their deployment, as well as to adverse environmental conditions [49]. Previously, a GWI animal model was established that attempted to capture the concomitant exposure to neurotoxic agents and high physiological stress, such as excessive heat and exercise, that was likely experienced during deployment by many of the veterans with GWI. By combining exposure to the sarin nerve agent surrogate disopropyl fluorophosphate (DFP) with prior subchronic exposure to the stress hormone corticosterone (CORT), it was shown that these exposures produce significantly exacerbated neuroinflammation and pain responses, as well as changes in the structure and cellular biology of the brain and glia [50-61]. While these results have been found acutely following exposure, they have begun to be supported in current evaluations of veterans with GWI by various imaging methods, finding evidence of both neuroinflammation [13] and brain structural changes indicative of the restructuring of the glia [62]. As such, it was important to investigate the impacts of these initial CORT and DFP exposures at longer-term time points that are more relevant to the current condition of veterans with GWI. Thus, to further represent the protracted nature of GWI (i.e., 30 years since the initiating event), mice were aged for a total of 11 weeks with either H 2 O or repeated exposure to 4-day bouts of CORT every other week to reflect ongoing stress, with the potential to be representative of PTSD, following the initial exposure to CORT and DFP as utilized in previous studies ( FIG. 2A ) [50-61]. Following this, a lipopolysaccharide (LPS) challenge was used to stimulate the immune system in lieu of exercise. To gauge the effect of the DFP exposure and LPS challenge, saline solution (SL) was used as a control. This resulted in five mouse groups analogous to the human cohort (H2O SL SL=HC; H 2 O DFP SL=GWI L  T0; CORT DFP SL=GWI H  T0; H 2 O DFP LPS=GWI L  T1; and CORT DFP LPS=GWI H  T1). 
     Serum cytokine levels for the GWI mouse model groups as compared to the H 2 O SL SL control are presented in  FIG. 2B . Statistical comparisons were made using heteroscedastic two-tailed t-tests between all GWI mouse model groups and the H 2 O SL SL control for both rest and LPS immune challenge timepoints, and within challenge groups for the H2O and CORT conditions. While the majority of measures for the unchallenged groups show no significant differences, there are statistical differences between the unchallenged groups and the control for IL-1α and IL-1β, and between the unchallenged groups for IL-4 and IL-5. 
     Turning to  FIG. 3A  through  FIG. 3C , shown is cytokine expression in a long-term mouse model of Gulf War Illness. (A) Graphical representation of the animal dosing paradigm. Mice were exposed to corticosterone (CORT; 200 mg/L in 0.6% ethanol) for 4 days at a time every other week for a total of 11 weeks. DFP (4 mg/kg, i.p.) was given on Day 5 and the LPS challenge (0.5 mg/kg, s.c.) was given after the final week of CORT. (B) Individual serum cytokine values normalized to the water (H2O), saline solution (SL), and SL condition values. Bars represent the mean ratios±SEM for each cytokine. Significance was set at p&lt;0.05 compared to H2O SL (*), and within challenge groups (&amp;: H2O DFP SL vs. H2O DFP LPS; and #: CORT DFP SL vs. CORT DFP LPS). The false discovery rate, as calculated via Storey et al. [48], was &lt;0.1 for all significant measures. (C) Cytokine mRNA was measured in the frontal cortex (n=2-5/group) by qRT-PCR 6 h following exposure to LPS. Bars represent the mean±SEM. Significance was set at p≤0.05 for comparisons within challenge group (*: DFP vs. CORT DFP) or across challenge group (#: SL vs. LPS). 
     The starkest differences are for the CORT DFP LPS group compared to the control and the non-stressed H2O DFP LPS group for IL-1α, IL-1β, IL-5, IL-6, KC (mouse equivalent of IL-8), mIL-10, IFNγ, and TNFα. While trends appear for mIL-6 and KC in the H2O DFP LPS group to be increased over the control, statistical difference is not reached due to large variance and a small group size. 
     Brain cytokine mRNA levels for the GWI mouse model groups as compared to the H 2 O SL control are presented in  FIG. 3C . Here, again, it is seen that the CORT DFP LPS group has statistically higher amounts of the proinflammatory cytokines IL-1β, IL-6, and TNFα□ than both the control and the other mouse groups. Additionally, this group also has elevated markers for microglial activation C—C Motif Chemokine Ligand 2 (CCL2) and oncostatin M (OSM). No other measures showed significant increases from the control or within groups. 
     3. Male BBB Neuro-Endocrine-Immune Model 
     Based on our previous work, we have combined a simple model of the overarching HPA-HPG-Immune system [38] with a more complex description of the immune function [39] and a brain compartment describing the neuroimmune interaction of the neurons and glia [53] ( FIG. 4 ). 
     Turning to  FIG. 4 , shown is a theoretical male neuroimmune signaling network. Light blue nodes denote the HPA axis model described in [37]. Dark blue nodes denote the male HPG axis as described in [38]. Orange and green nodes on the right denote a system of immune signaling molecules and immune cells, respectively, as originally described in [39]. Grey, yellow, pink, purple, and green nodes on the left denote neuronal-glia interactions, growth factors, blood-brain barrier, neurotransmitter, and neuroimmune signaling molecules, respectively, as originally described in [53]. Red nodes denote external influences on the system. Green edges are stimulatory, and red edges are inhibitory. 
       FIG. 5  illustrates predicted stable states of the theoretical male neuroimmune signaling network presented in  FIG. 4 . White fill indicates a nominal state the same as the healthy control (0); green a state higher than the healthy control (1); and red state lower than healthy control (−1). 
     Using ternary logic modeling and performing Monte Carlo simulations of state evolution in this theoretical male neuroimmune signaling network identified six potential stable behaviors available to the system ( FIG. 5 ). The reference stable state (SS0) is that of typical health in which all measures are nominally the same as the healthy control condition. The first alternate stable state (SS1) displays a low adrenocorticotropic hormone (ACTH) concentration with a high glucocorticoid receptor number in both its single membrane bound (GR) and internalized dimerized (GRD) forms. Although Monte Carlo algorithms are used, it is important to note that other computational algorithms such as a genetic algorithm, an optimization algorithm, and evolutionary algorithm, or a combination thereof with Monte Carlo are possible. 
     The remaining alternate states (SS2, SS3, SS4, and SS5) all present with elevated HPA axis activity: cortisol, GR, and GRD, and depressed HPG activity: testosterone (TEST), luteinizing hormone (LH)/follicle-stimulating hormone (FSH), gonadotropin releasing hormone (GnRH) in addition to depressed natural killer cell activity (NK), and an increase in monokine (MK) 2 activity (IL-10, IL-22). Beyond this, the second alternate state (SS2) is distinguished by a periphery marked with elevated chemokines (CK) 2 (interleukin (IL)-4, IL-5, IL-13), CK17 (IL-17A, IL-17F), MK1 (IL-1α, IL-1β, IL-8, IL-12), MK6 (IL-6), MK 21 (IL-21), T-helper (Th) 2 cells, Th17B and Th17(23), and depressed MK15 (IL-15), MK23 (IL-23), MK27 (IL-27), dendritic cells (DC), and T regulatory (Treg) cells. In the brain compartment, SS2 presents elevated endothelial cell activity, increased T-cell activity, IL-4, and decreased pro-inflammatory cytokines (IL-1β, IL-6, TNFα), and insulin-like growth factor 1 (IGF-1). The third alternate state (SS3) presents in the periphery with elevated CK1 (IL-2, IFNγ, TNFα, TNFβ), cytotoxic lymphocytes (CTL), MK2, and Tregs, and in the brain with elevated endothelial cell activity, pro-inflammatory cytokines, and T cell activity, and depressed IL-4 and IGF-1. The fourth alternate state (SS4) likewise has elevated CK1, CTL, MK2 endothelial cell activity in the periphery, and elevated IL-4 and T-Cells in the brain. Additionally, SS4 displays elevated CK17, MK1, MK6, MK21, Th17B, Th17(23), and depressed transforming growth factor β (TGFβ) and Th1 in the periphery. The final alternate state (SS5) is marked by a periphery with elevated CK1, CK17, MK1, MK2, MK6, MK15, CTL, Th17B, and Th17(23), and in the brain with elevated pro-inflammatory cytokines, microglia activity, and vascular endothelial growth factor (VEGF), and depressed endothelial cell function. 
     4. Alignment of the Blood Profiles with the Model-Predicted States 
     To determine which of the animal model conditions best represents the GWI cohort subgroups, to ultimately infer the brain state the GWI cohort subgroups, the above model describing the known interactions of the neuroimmune system was used to provide a context for comparisons between human and animal data. To do this, each subgroup condition in both the human and animal groups was compared to their respective control condition to determine the degree of agreement with each stable behavior predicted by the model ( FIG. 5 ). Initially, blood profiles for the baseline clinical data obtained from healthy control patients and data for aged mice receiving a course of pure drinking water, an initial saline injection, and a follow-up injection of saline were used to define the base topology of the simulated stable states using Brown&#39;s method, as described in the Materials and Methods section ( FIG. 6A  through  FIG. 6D ). As can be seen from  FIG. 6A  and  FIG. 6B , the overall topology of the stable state locations is retained between humans and mice, as evidenced by the relative spacing between points. However, there is a discrepancy in the actual distance values between the stable points, owing to the difference in the human and mouse systems. By comparing the edge-by-edge distances between these networks, a strong linear correlation between mouse and human data is revealed, which is used to scale the mouse net-work data to humans for further analysis ( FIG. 6C  and  FIG. 6D ). 
     More specifically,  FIG. 6  is a reference blood profile topological network. (A) Topology of the distance between the simulated stable states using Brown&#39;s method based on human healthy control data. (B) Topology of the distance between the simulated stable states using Brown&#39;s method based on data from mice receiving a course of pure drinking water, an initial saline injection, and a follow-up injection of saline. (C) Regression of the human topology onto the mouse topology. (D) Topology of the distance between the simulated stable states using Brown&#39;s method based on data from mice receiving a course of pure drinking water, an initial saline injection, and a follow-up injection of saline scaled to the human topology. 
     Turning now to  FIG. 7 , shown is a comparison of the GWI trauma subtype blood profile topologies with mouse models at rest. (A) Topology of the distance between the simulated stable states for GWI L  at rest. (B) Topology of the distance between the simulated stable states for GWI H  at rest. (C) Topology of the distance between the simulated stable states from mice receiving a course of pure drinking water, an injection of DFP, and a follow-up injection of saline. (D) Topology of the distance between the simulated stable states from mice receiving a course of drinking water with CORT, an injection of DFP, and a follow-up injection of saline. 
     Cohort and mouse model subgroup condition blood profiles at rest were then added to the reference blood profile topological networks ( FIG. 7A  through  FIG. 7D ). The topological profile of the resulting network for GWI H  ( FIG. 7B ) differed from GWI L  ( FIG. 7A ) with an RMSD of 0.42. Similar topological differences are seen for DFP-exposed mice ( FIG. 7D ) with and without prior treatment with CORT ( FIG. 7C ). The difference between these mice conditions is comparable to the human conditions given by an RMSD of 0.33. 
     Comparison of the GWI illness conditions to the mouse conditions at rest show that the GWI L  condition is most closely matched by the non-CORT-treated mouse model (RMSD=0.06) rather than the CORT-treated model (RMSD=0.36). The opposite was shown for the GWI H  condition, with the CORT-treated model (RMSD=0.28) being a better match over the non-CORT-treated model (RMSD=0.47). However, despite this agreement between the trauma/stressed states, there are noticeable differences that account for the non-zero RMSD. Both GWI illness subtypes show the closest alignment with stable state SS4/SS5. This difference is driven by GWI H  having a large distance from SS0/SS1 and closer alignment to SS3, while this is reversed for GWI L . 
     This is not the case for the mouse models as the CORT-treated model shows the closest alignment with SS3, while the non-CORT-treated model shows the closest alignment with SS0/SS1. Yet, while this is a noticeable difference, the relationship between the illness condition and all four stable states determines its overall placement in the landscape and the resulting topology. 
     Turning now to  FIG. 8  is a comparison of the GWI trauma subtype blood profile topologies with mouse models under challenge. (A) Topology of the distance between the simulated stable states for GWIL at peak exercise. (B) Topology of the distance between the simulated stable states for GWIH at peak exercise. (C) Topology of the distance between the simulated stable states of mice receiving a course of pure drinking water, an injection of DFP, and a follow-up injection of LPS. (D) Topology of the distance between the simulated stable states of mice receiving a course of drinking water with CORT, an injection of DFP, and a follow-up injection of LPS. 
     Blood profiles for the model and illness conditions under challenge were separately added to the healthy blood profile topological networks ( FIG. 8 ). Unlike the at-rest condition, the non-trauma GWIL and trauma GWIH under challenge conditions ( FIG. 8A  and  FIG. 8B , respectively) present with very similar topologies (RMSD=0.02), as do the non-CORT- and CORT-treated DFP LPS mouse models, although to a lesser degree (RMSD=0.05). Comparison of GWI illness conditions to the mouse conditions under challenge shows that the GWIL condition is still most closely matched by the non-CORT treated mouse model (RMSD=0.02) ( FIG. 8C ) rather than the CORT-treated model (RMSD=0.03) ( FIG. 8D ). However, the degree of difference is much less. Again, the GWIH condition showed a better match with the CORT-treated model (RMSD=0.02), being a better match over the non-CORT-treated model (RMSD=0.04), again with smaller differences. Unlike the at-rest conditions, the agreements between the trauma/stressed states share the same ordering of differences between the stable state positions with only minor differences to account for the small non-zero RMSDs. In all cases, the conditions show the closest alignment with stable state SS4/SS5, followed by SS2, SS3, and finally SS0/SS1. 
       FIG. 9A  through  FIG. 9E  are a comparison of the mouse model conditions brain and blood profile topologies. (A) Topology of the distance between the simulated stable states using Brown&#39;s method based on data from mice receiving a course of pure drinking water, an initial saline injection, and a follow-up injection of saline. (B) Topology of the distance between the simulated stable states from mice receiving a course of pure drinking water, an injection of DFP, and a follow-up injection of saline. (C) Topology of the distance between the simulated stable states from mice receiving a course of drinking water with CORT, an injection of DFP, and a follow-up injection of saline. (D) Topology of the distance between the simulated stable states of mice receiving a course of pure drinking water, an injection of DFP, and a follow-up injection of LPS. (E) Topology of the distance between the simulated stable states of mice receiving a course of drinking water with CORT, an injection of DFP, and a follow-up injection of LPS. 
     As it is not possible to assess brain measures in the human cohort directly, adding brain measures from the mouse models to their blood profiles allows for the discrimination between SS4 and SS5 while retaining the overall structure observed in the reference brain and blood profiles ( FIG. 9A ), and therefore infer the probable brain states for the human cohort groups. Topologically, the non-challenged mouse model, stressed and un-stressed, resembles the blood-only profiles, as do the challenged stressed and unstressed conditions. For the non-challenged conditions ( FIG. 9B ,  FIG. 9B  C), the Water DFP SL mouse model retains its closest proximity to SS0/SS1, while the CORT DFP SL mouse model retains closest proximity to SS3. For the challenged conditions ( FIG. 9D ,  FIG. 9E ), the Water DFP LPS mouse model is most comparable to the SS5 state, while the CORT DFP LPS mouse model is equidistant from both SS4 and SS5, with only a slightly longer distance to SS2. The remainder of the edge orders for all models is also retained from the blood-only comparisons ( FIG. 7C ,  FIG. 7D  and  FIG. 8C ,  FIG. 8D ), save for SS5. 
     5. Simulated Treatment Courses 
     As determined above, the GWI blood profiles both align more closely with the SS4/SS5 attractor; however, the blood profiles alone cannot discriminate between these two attractor states. The alignment of the mouse profiles, GWIL, and GWIH are best represented by the Water-DFP and CORT-DFP models, respectively, both at rest (saline) and under LPS challenge. The brain profiles of these mouse models suggest that GWIL orbits SS5, while GWIH orbits equidistant from states SS4 and SS5 closely followed by SS2, all of which are taken as starting points for the simulation of the treatments. Treatments were simulated, allowing for increases in cortisol, ACTH, and testosterone, inhibition of glucocorticoid receptors, and bidirectional modulation of all peripheral cytokines. 
     Finally, for the SS2 start state, only two interventions were needed to obtain at least a 1% return to health ( FIG. 12  with Table 4). This two-intervention treatment course consisted of a near overlapping continual treatment of inhibiting pro-inflammatory cytokine IL-6 with inhibition of the glucocorticoid receptor. This dual combination strategy resulted in a 31% return to health, less than both the SS4 and SS5 three-target treatments, but with less intervention overall. 
     For the SS5 start state, a minimum of three treatment targets were needed to obtain at least a 1% return to health ( FIG. 10  with Table 2). The treatment course consisted of an initial phase of increasing testosterone, followed by an inhibition of CK1 cytokines, including IL-2, IFNγ, TNFα, and TNFβ, followed by inhibiting the glucocorticoid receptor with some overlap between these phases. This minimum intervention resulted in a 37% return to health. 
     For the SS4 start state, again, a minimum of three treatment targets were needed to obtain at least 1% return to health ( FIG. 11  with Table 3). The treatment course consisted of a continual treatment of inhibiting the pro-inflammatory cytokine IL-6, which contained an initial phase of inhibiting CK1 cytokines including IL-2, IFNγ, TNFα, and TNFβ, followed by a phase of inhibiting the glucocorticoid receptor with no overlap between these two contained phases. This minimum intervention resulted in an improved 61% return to health over the SS4 prescribed treatment course. 
     Discussion 
     The wide range of symptoms observed in GWI is consistent with underlying neuroinflammatory processes. Ultimately, these processes have effects beyond the CNS, affecting immune and endocrine function in the periphery and ultimately giving rise to the symptom presentations in the illness. As the presence of PTSD symptoms exacerbates the overall clinical symptoms of GWI [19], and PTSD shows its immunological dysfunction [21-24], it is reasonable to assume that past trauma may be used to delineate two biological GWI subgroups. In the context of the neuroimmune model presented here, the high and low trauma GWI group blood profiles at rest share the closest proximity with the same altered regulatory state marked by elevated HPA axis activity, depressed HPG activity, depressed NK activity, and increased MK2, CK1, CK17, MK1, MK2, MK6, MK15, CTL, Th17B, and Th17(23) activity, albeit not with complete agreement. Under exercise challenge, these two groups converge to this common point and are nearly identical. 
     These results alone suggest that the absence and presence of PTSD do not delineate specific biological subtypes of GWI. However, when projected into the neuroimmune model state space, these two subtypes present as topologically distinct at rest, suggesting that the underlying biological regulation between these two groups are unique. These modeling results are consistent with past research into the biomarkers/biosignatures of GWI. The difference in biological profiles between GWI and the controls has been shown to be more distinct and consistent under an exercise challenge [6]. At rest, GWI profiles are variable and sometimes present with inconsistent results [63,64]. Past biomarker/biosignature studies have not grouped GWI in terms of prior trauma; as such, the distinct GWI L  and GWI H  subgroups intermixed at rest would be expected to yield mixed and inconsistent results compared to controls as the proportion of each group varies. However, under exercise challenge, as the profiles of GWI L  and GWI H  converge, the typical neuro-endocrine immune signature emerges as different from the control. 
     The difference between the profiles of GWI L  and GWI H  at rest is due to their relation to the remaining stable regulatory patterns that define their difference in blood alone. As neither of the subgroups aligns completely with any of the defined stable states, there are factors not accounted for in the model to hold the subgroups at their baseline positions. This may be in both the periphery or in the brain portions of the network model. Further refinement of the granularity of the peripheral endocrine and immune system elements is thus hypothesized to further define these two subgroups in human subjects. 
     The topological alignment between the GWI subtypes and the animal models suggests that the use of corticosterone priming of the neuroinflammatory response seems to best represent GWI with past traumatic events (i.e., GWI with probable PTSD; GWI H ), while the animal model without corticosterone priming best represents GWI alone (i.e., GWI without probable PTSD; GWI L ). However, while the animal models at rest topologically resemble the GWI H  and GWI L  states, the challenged condition in human and nonhuman animals is more closely aligned and represents the best scenario for assessing treatments. 
     By applying our genetic algorithm methodology to these starting conditions, we identify three possible candidate treatment courses for GWI. The first, involving increasing testosterone, followed by inhibiting IL-2, IFNγ, TNFα, or TNFβ, followed by inhibiting the glucocorticoid receptor, applies to both the GWI L  and GWI H  subgroups and is the only predicted treatment course for GWI L  based on this model. This treatment course is consistent with past findings regarding Th1 immunity and glucocorticoids in GWI and our previously predicted treatment course for GWI involving modulating Th1 immunity via inhibiting TNFα and TNFβ via etanercept, followed by blocking the glucocorticoid receptor with mifepristone [65,66], except for the addition of the testosterone pre-treatment. 
     There are many anecdotal reports that male veterans with GWI suffer from low levels of testosterone; however, there is a lack of formal studies investigating this claim. Studies of self-reported sexual dysfunction have included reports of decreased libido, erectile dysfunction, discomfort or pain during intercourse, and a burning sensation after sex. A primary study on sexual dysfunction in a Danish cohort of Gulf War veterans found that these sexual problems were associated with traumatic events and not related to differences in major reproductive hormone parameters [67,68]. However, these studies concern Gulf War veterans, which who do not necessarily meet the criteria for Gulf War Illness. Furthermore, Danish troops were in the Persian Gulf region after the war as peacekeepers and were not exposed to sarin [69]. Other studies investigating the testosterone levels in Gulf War veterans compare groups with high or low exposure to depleted uranium. These studies have similarly found no statistically significant differences between these groups in endocrine function (i.e., including serum follicle-stimulating hormone, luteinizing hormone, prolactin, and total testosterone) [70-76]. This, however, does not rule out the possibility of low testosterone in veterans with Gulf War Illness compared to sedentary Gulf War-era veteran controls. Several studies in rats have identified anti-androgenic effects (i.e., reduced levels of testosterone, follicle-stimulating hormone, and luteinizing hormone, as well as testicular shrinkage) of the GW-relevant organophosphate acetylcholinesterase inhibitor and pesticide, chlorpyrifos [77-81]. Studies investigating the diffuse chronic pain of fibromyalgia, which has commonalities with other chronic pain conditions such as PTSD and Gulf War Illness, have suggested low or deficient testosterone serum levels that likely cause an inflamed nociceptive nervous system [82]. This hypothesis is supported by the relief of chronic pain symptoms following the application of testosterone gel [83]. Therefore, while our modeling efforts indicate the potential for low testosterone in GWI [38,39,66,84], which is supported by trends in our data, there is a lack of conclusive findings in this area; thus, a more formal analysis of endocrine dysfunction in GWI accounting for confounding factors, such as age and BMI, is needed. 
     Beyond the first three-intervention treatment course predicted for GWI in total, our modeling efforts also identified two additional treatment courses for GWI H , each with its separate advantages and disadvantages. The second course is another three-intervention treatment strategy that involves a continuous inhibition of IL-6, which contains an inhibition of IL-2, IFNγ, TNFα, or TNFβ, followed by inhibition of the glucocorticoid receptor with no overlap. Like the addition of testosterone suggested for the first treatment course, this treatment course is consistent with the previously reported treatment course of modulation of Th1 immunity followed by blocking the glucocorticoid receptor, except here it involves the addition of a continual inhibition of IL-6 during the treatment course. Here, the potential drawback is both the continual administration of a potent IL-6 immune inhibitor over the treatment course in addition to another interspersed IL-2, IFNγ, TNFα, or TNFβ immune inhibitor. The advantage of this three-intervention treatment is the suggested doubling of the overall return to health rate compared to the first. 
     Finally, the third treatment course is a two-intervention strategy involving simultaneous inhibition of IL-6 and the glucocorticoid receptor. This is a novel treatment course for GWI H  of minimal intervention with a comparable return to the health rate of the first treatment course. It does involve the continual administration of a potent IL-6 immune inhibitor over the treatment course in conjunction with glucocorticoid inhibition but does not have the same drawback as the second treatment course involving simultaneous dual immune inhibition. 
     The inhibition of IL-6 for the treatment of GWI with probable PTSD makes sense in the light of the increased levels of IL-6 found in PTSD. For example, elevated serum IL-6 and IL-6 receptor concentrations have been found in PTSD following accidental man-made traumatic events [85], and enhanced IL-6 response to the mental stress following myocardial infarction [86]. A meta-analysis and meta-regression of 8057 studies comparing inflammatory markers between patients with PTSD and healthy controls found increased levels of IL-6 that are more pronounced in PTSD without major depressive disorder and for those not taking medication, with the severity of PTSD being positively associated with IL-6 levels [87]. Many animal models have also been used to understand how the brain responds to trauma and investigate treatments that can reduce symptoms of PTSD. Despite the varied methodologies used to model PTSD, rodent models have consistently shown elevated inflammatory markers, including IL-6, both in the central nervous system and peripherally after exposure to inescapable shock [88], chronic variable stress [89,90], and social defeat [91]. This indicates that IL-6 is a feature of PTSD-related pathology both in human patients and in rodent models. Interestingly, treatments that alleviate PTSD-related symptoms also reduce IL-6. This finding is corroborated with anti-inflammatory drugs such as minocycline, which reduces IL-6 and anxiety-related behaviors in rodent models of PTSD [92,93]. Not only do these anti-inflammatory drugs have anxiolytic effects, indicating that reducing pro-inflammatory cytokines could ameliorate PTSD, but also the frontline drugs used to treat PTSD, serotonin selective reuptake inhibitors (SSRIs), reduce inflammatory markers such as IL-6 [92,94]. Taken together, a trauma that induces PTSD-like features in nonhuman animals increase IL-6 and treatments that reduce anxiety after trauma also reduce IL-6, so it follows that IL-6 may play a role in the mechanism of PTSD, but most data supporting this role of IL-6 are indirect. More recently, inhibiting IL-6 using antibodies or transgenic silencing showed that IL-6 plays a causal role in the enhanced fear memory commonly seen in PTSD [95]. Thus, animal models of PTSD indicate that inflammatory markers, including IL-6, are both a feature and a mechanism underlying PTSD. 
     Materials and Methods 
     Discrete Ternary Logical Analysis 
     The discrete ternary logical network analysis used in the present work is an extension of a methodology proposed by Mendoza and Xenarios [40] and Thomas [43], and has been reported previously by our group [38,39,53,66]. We encode documented feedback mechanisms within the endocrine-immune system using only the direction (source and target) and type (activator or inhibitor) of interaction. As data describing the magnitude of changes remains limited, we consider all cell types to be equally responsive to the actions of the cytokines for which they express receptors. Accordingly, we also consider cytokine synthesis to be equivalent regardless of cell type. Using this formalism, we determine the number and type of stable resting states supported by the regulatory circuitry as well as the specific qualitative endocrine-immune signatures at each of these stable points without requiring detailed kinetic information; that is, we determine where the system would eventually come to rest even though we may not know how quickly this equilibrium will be reached. 
     In this model, signaling molecules and cell types are represented as individual variables, each capable of adopting three discrete states: −1 (downregulated), 0 (nominal), and 1 (upregulated). At any point in time t, the state of a system with n variables can be represented by the vector {right arrow over (x)}(t), such that 
         {right arrow over (x)} ( t )=( x   1 ( t ), x   2 ( t ), . . . , x   N ( t )),  (1)
 
     where x i (t) is the state of the i th  variable of the n variable system at time t. The image vector {right arrow over (x)}(t)(t+1) describes the preferred state towards 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 a variable and the mode of action (i.e., activate or inhibit) of the neighboring input variables. These logic statements are expressed as follows (Equation (2)): 
     
       
         
           
             
               
                 
                   
                     
                       x 
                       i 
                     
                     ⁡ 
                     
                       ( 
                       
                         t 
                         + 
                         1 
                       
                       ) 
                     
                   
                   = 
                   
                     { 
                     
                       
                         
                           
                             
                               
                                 ( 
                                 
                                   
                                     
                                       x 
                                       
                                         i 
                                         ⁢ 
                                         1 
                                       
                                       A 
                                     
                                     ⁡ 
                                     
                                       ( 
                                       t 
                                       ) 
                                     
                                   
                                   ⋁ 
                                   
                                     
                                       x 
                                       
                                         i 
                                         ⁢ 
                                         2 
                                       
                                       A 
                                     
                                     ⁡ 
                                     
                                       ( 
                                       t 
                                       ) 
                                     
                                   
                                   ⋁ 
                                   
                                       
                                   
                                   ⁢ 
                                   … 
                                   ⁢ 
                                   
                                       
                                   
                                   ⋁ 
                                   
                                     
                                       x 
                                       
                                         i 
                                         ⁢ 
                                         X 
                                       
                                       A 
                                     
                                     ⁡ 
                                     
                                       ( 
                                       t 
                                       ) 
                                     
                                   
                                 
                                 ) 
                               
                               ⁢ 
                               
                                 ∇ 
                                 
                                   ( 
                                   
                                     
                                       
                                         x 
                                         
                                           i 
                                           ⁢ 
                                           1 
                                         
                                         I 
                                       
                                       ⁡ 
                                       
                                         ( 
                                         t 
                                         ) 
                                       
                                     
                                     ⋁ 
                                     
                                       
                                         x 
                                         
                                           i 
                                           ⁢ 
                                           2 
                                         
                                         I 
                                       
                                       ⁡ 
                                       
                                         ( 
                                         t 
                                         ) 
                                       
                                     
                                     ⋁ 
                                     
                                         
                                     
                                     ⁢ 
                                     … 
                                     ⁢ 
                                     
                                         
                                     
                                     ⋁ 
                                     
                                       
                                         x 
                                         
                                           i 
                                           ⁢ 
                                           Y 
                                         
                                         I 
                                       
                                       ⁡ 
                                       
                                         ( 
                                         t 
                                         ) 
                                       
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                         
                         
                           
                             
                               ( 
                               
                                 
                                   
                                     x 
                                     
                                       i 
                                       ⁢ 
                                       1 
                                     
                                     A 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     t 
                                     ) 
                                   
                                 
                                 ⋁ 
                                 
                                   
                                     x 
                                     
                                       i 
                                       ⁢ 
                                       2 
                                     
                                     A 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     t 
                                     ) 
                                   
                                 
                                 ⋁ 
                                 
                                     
                                 
                                 ⁢ 
                                 … 
                                 ⁢ 
                                 
                                     
                                 
                                 ⋁ 
                                 
                                   
                                     x 
                                     
                                       i 
                                       ⁢ 
                                       X 
                                     
                                     A 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     t 
                                     ) 
                                   
                                 
                               
                               ) 
                             
                           
                         
                         
                           
                             
                               ⫬ 
                               
                                 ( 
                                 
                                   
                                     
                                       x 
                                       
                                         i 
                                         ⁢ 
                                         1 
                                       
                                       I 
                                     
                                     ⁡ 
                                     
                                       ( 
                                       t 
                                       ) 
                                     
                                   
                                   ⋁ 
                                   
                                     
                                       x 
                                       
                                         i 
                                         ⁢ 
                                         2 
                                       
                                       I 
                                     
                                     ⁡ 
                                     
                                       ( 
                                       t 
                                       ) 
                                     
                                   
                                   ⋁ 
                                   
                                       
                                   
                                   ⁢ 
                                   … 
                                   ⁢ 
                                   
                                       
                                   
                                   ⋁ 
                                   
                                     
                                       x 
                                       
                                         i 
                                         ⁢ 
                                         Y 
                                       
                                       I 
                                     
                                     ⁡ 
                                     
                                       ( 
                                       t 
                                       ) 
                                     
                                   
                                 
                                 ) 
                               
                             
                           
                         
                       
                       , 
                     
                   
                 
               
               
                 
                   ( 
                   2 
                   ) 
                 
               
             
           
         
       
     
     where the ∇, v, and ¬ symbols are ternary HIGH/LOW PASS, OR and NOT operators, x ij   A  is the state of the i th  variable&#39;s j th  activator, and x ij   I  is the state of the i th  variable&#39;s k th  inhibitor. The ternary operators given in Equation (2) are described in further detail in [38]. 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. While the number of X activators and/or Y inhibitors for a given variable may remain static, they may also be allowed to change based on predefined conditions, such as the state of one or more variables. In the case where the number of X activators and/or Y inhibitors are conditionally dependent, Equation (2) is modified to reflect the number of activators or inhibitors allowed by the predefined conditions and the current state of the network system. In this case of this model conditional edges are only allowed if the Endothelial Cells node is −1. 
     2. Monte Carlo Simulation of State Evolution 
     Following [38], the evolution of state transitions supported by the model was analyzed by developing a Monte Carlo simulation algorithm. From any initial starting state, allowable state transitions are determined based on Equation (2). Applying Equation (2) to each variable in the model for the state of the system at time t, {right arrow over (x)}(t) defines the image vector {right arrow over (x)}(t+1) for the next timepoint. With {right arrow over (x)}(t+1) defined, the system is updated asynchronously allowing only one variable to change at each timestep following the generalized logical analysis of Thomas [43]. According to this method, one variable slated to change by the image vector is chosen at random using a uniform equal distribution and used to generate the next allowable state; that is, if the i th  variable of the state vector, x i (t), is chosen at random to evolve it is moved one step towards its preferred image x i (t+1) (e.g., if x i (t)=−1 and x i (t+1)=1, then x i (t+1) is set to 0). Thus, for each current state of the system there are potentially several subsequent states towards which it may asynchronously evolve. States for which the image vector is the same as the current state vector are considered stable (steady states, attractors, basins, etc.), and do not evolve further in time. The Monte Carlo procedure is performed until such a stable state is reached. Executing the simulation multiple times gives a distribution of paths that is used to determine the behavior of the system from any given start state. 
     3. Participants/Procedure 
     All participants signed an informed consent approved by the Institutional Review Board of the Miami Veterans Affairs Medical Center. Ethics review and approval for data analysis was also obtained by the IRB of Nova Southeastern University. 
     Participants were recruited via the Miami Veterans Affairs (VA) Medical Health Center in two cohorts funded under a Veterans Affairs Merit award (GWI: n=27, healthy controls (HC): n=19) and a Department of Defense Gulf War Illness Research Program award (W81XWH-09-2-0071) (GWI: n=24, HC: n=19), both of which compared male veterans with GWI to HC. Therefore, the full sample for the multivariate and univariate analyses included n=51 male veterans with GWI in addition to n=38 healthy controls. 
     Inclusion criteria for GWI participants was derived from Fukuda et al. [96] and consisted in identifying veterans deployed to the theater of operations between 8 Aug. 1990, and 31 Jul. 1991, with one or more symptoms present for 6 months from at least 2 of the following: fatigue, mood, and cognitive complaints, as well as musculoskeletal complaints. Participants were in good health prior to 1990 and had no current exclusionary diagnoses defined by Reeves et al. [97]. This includes the exclusion of major dementias of any type and alcoholism or drug abuse, medical conditions including organ failure, rheumatologic disorders, and use of medications that impact immune function, such as steroids or immunosuppressants. Collins et al. [98] supports the use of the Fukuda definition in GWI. Control participants consisted of Gulf War era veterans self-defined as healthy with no exclusionary diagnoses and sedentary (no regular exercise program, sedentary employment). 
     All participants were subjected to a standard maximal graded exercise test to stimulate their immune response, and blood samples were collected at rest (T0) and peak effort (T1) as previously described [6]. 
     4. Human Blood Analysis 
     Plasma was collected immediately upon receiving. EDTA Plasma was collected by centrifugation, 2000×g for 15 min, and frozen immediately. Cytokine analysis was performed using Quansys chemiluminescent assays (Quansys Biosciences, Logan, Utah). The Quansys Imager, driven by an 8.4-megapixel Canon 20D digital SLR camera, supports 96-well plate-based chemiluminescent imaging. The Q-Plex™ Human Cytokine Screen is a quantitative ELISA-based test where distinct capture antibodies have been absorbed to each well of a 96-well plate in a defined array. Custom 16- and 18-multiplex assays were used in the Veterans Affairs Merit award and a Department of Defense Gulf War Illness Research Program award (W81XWH-09-2-0071) studies, respectively. 
     The 16-multiplex panel includes TNFα, TNFβ, IL-1α, IL-1β, IL-6, IFN-γ, IL-12, IL-2, IL-15, IL-8, IL-5, IL-17, IL-23, IL-10, and IL-13 (Fletcher et al. (2009)). For the standard curves, we used the second-order (k=2) polynomial regression model (parabolic curve), Y=b 0 +b 1 X+b 2 × 2  . . . +b k X k , where Y caret is the predicted outcome value for the polynomial model with regression coefficients b 1  to k for each degree and y intercept b 0 . Quadruplicate determinations were made, i.e., each sample was run in duplicate in two separate assays. 
     The 18-multiplex panel includes TNFα, TNFβ, IL-1α, IL-1β, IL-6, TNF-RI, TNF-RII, IFN-γ, IL-12, IL-2, IL-15, IL-8, IL-5, IL-17, IL-23, IL-10, and IL-13. Briefly, plasma samples are thawed at 4° C. overnight. Samples are plated in duplicate following the manufacturer&#39;s protocol. The plates were read at 270 s using the Q-view Imager LS (Quansys). Individual cytokine concentrations were obtained using image analysis software (Q-view, Quansys Biosciences). Sample concentrations were calculated from standard curves created by a five-parameter logistic regression (5PL) with √y weighting. The average value from each duplicate was then used for subsequent analyses. 
     To account for 16-plex and 18-plex differences, data were separated and standardized by z-scoring. The two data sets were then combined and reverse z-scored using the mean and standard deviation of the 18-plex dataset to scale the 16-plex data to the 18-plex set. Following this, the dataset was normalized using min-max normalization so the final data values ranged between 0 and 1. Normalized cytokine values are given in Table 1 of  FIG. 1 . 
     5. Animals and Dosing 
     All animal experiments were performed using protocols approved by the Centers for Disease Control—National Institute for Occupational Safety and Health Institutional Animal Care and Use Committee (CDC-Morgantown IACUC) and the United States Army Medical Research and Development Command Animal Care and Use Review Office (USAMRDC ACURO). Adult male C57BL/6J mice (8-12 weeks of age; weighing approximately 22 g) were purchased from Jackson Labs (Bar Harbor, Me., USA). Mice were single housed in a temperature- (21±1° C.) and humidity-controlled (50±10%) room under filtered positive-pressure filtration, and a 12 h light/dark cycle (lights on at 0600 EDT); mice were given ad libitum access to food (Harlan 7913 irradiated NIH-31 modified 6% rodent chow) and water. 
     For dosing, mice were randomly placed into groups (n=5/group) and either remained on regular drinking water or pretreated with CORT (Steraloids Inc., Newport, R.I., USA) in the drinking water (200 mg/L in 0.6% Ethanol (EtOH; Sigma-Aldrich, Inc., St. Louis, Mo., USA)) ad libitum for 4 days followed by exposure to DFP (4 mg/kg, i.p.; Sigma-Aldrich, Inc., St. Louis, Mo., USA) or saline (0.9%) on the 5th day. Following this exposure combination, modeling the initiating factors for GWI [50,51,54,61], mice were repeatedly exposed to 4-day bouts of CORT every other week (4 days on/10 days off) for a total of 11 weeks. Following the final 4-day CORT exposure, mice were challenged with the inflammatory bacterial endotoxin lipopolysaccharide (LPS; 0.5 mg/kg, s.c.; Sigma-Aldrich, Inc., St. Louis, Mo., USA) or saline and sacrificed by decapitation at 6 h following this exposure. Frontal cortex and blood serum samples were collected as previously described [50,51,54,61]. Due to mortality associated with DFP exposure, the DFP and DFP LPS groups comprised an n of 2 and 3, respectively. A schematic of the dosing paradigm is shown in  FIG. 3A . 
     6. Mouse Serum Cytokine Analysis 
     Cytokine levels were measured in mouse serum (n=2-5/group) as previously described [61]. Briefly, serum samples were analyzed using Q-Plex™ Mouse Cytokine Screen (16-plex), Quansys Imager, and Quansys reagents (Quansys Biosciences, Logan, Utah). Cytokine levels are expressed as pg/mL and presented as the mean±SEM in Supplementary Table S2. 
     7. qRT-PCR Analysis of Mouse Brain Cytokines 
     mRNA was isolated and cytokine expression measured by qRT-PCR analysis as previously described [50,99]. Briefly, total RNA was isolated from the frontal cortex using Trizol Reagent (Thermo Fisher Scientific, Waltham, Mass., USA), Phase-lock heavy gel (Eppendorf AG, Hamburg, Germany), and RNeasy mini spin columns (Qiagen, Valencia, Calif., USA) following the manufacturer&#39;s instructions. PCR analysis of the housekeeping gene, glyceraldehyde-3-phosphate dehydrogenase (GAPDH), and of the proinflammatory mediators, TNFα, IL-6, CCL2, IL-1β, leukemia inhibitor factor (LIF), and oncostatin M (OSM) was performed in an ABI7500 Real-Time PCR System (Thermo Fisher Scientific) in combination with TaqMan® chemistry. Relative quantification of gene expression was performed using the comparative threshold (ΔΔC T ) method to normalize expression changes against the GAPDH control, as well as normalize the expression changes to the corresponding saline-treated controls ( FIG. 3C ). Statistical significance (p&lt;0.05) was determined by two-way ANOVA followed by multiple comparisons using Fisher&#39;s Least Significant Difference in SigmaPlot (v. 14, Systat Software, Inc., San Jose, Calif., USA). 
     8. Comparison of the Gene Expression Data to the Model-Predicted States 
     Following [38], to check the validity of our model, we compared stable states predicted by the latter to gene expression profiles measured experimentally in the mouse model of Gulf War Illness [51] in both brain and blood, and to blood measures using the cohort of male veterans. Brown&#39;s theoretical approximation [100] of Fisher&#39;s statistics was used as a measure of similarity between a given model-predicted state and the expression profile measured in a subset of genes corresponding to the cellular and molecular entities used in our model as done in our previous work [38,39,66,84]. Fisher&#39;s statistics provide a meta-analysis technique to combine probabilities and obtain the overall significance of a set of p-values corresponding to independent tests of the same null hypothesis. The combined χ 2  statistic, 
     
       
         
           
             
               
                 
                   
                     T 
                     0 
                   
                   = 
                   
                     
                       - 
                       2 
                     
                     ⁢ 
                     
                       
                         ∑ 
                         
                           i 
                           = 
                           1 
                         
                         n 
                       
                       ⁢ 
                       
                         ln 
                         ⁡ 
                         
                           ( 
                           
                             p 
                             i 
                           
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
     where n is the number of measurable variables and p i  are the corresponding p-values under the null hypothesis, has a χ 2  distribution with 2N degrees of freedom, assuming the hypothesis tests are independent. As evidenced by the connectivity of the system studied here, these model entities do not express independently. As a result, the direct application of this test statistic is not valid since the assumption of independence is violated. To accommodate this, Brown [100] suggested a method for combining non-independent tests. If the tests are not independent, then the statistic T 0  has mean m=2N and variance (σ 2 ), given as 
       σ 2 =4 N+ 2Σ i=1   N−1 Σ j=i+1   N cov(−2 lnp   i ,−2 lnp   j )  (4)
 
     where p i  and p j  are the p-values for each test and the covariance (coy) is calculated as 
     
       
         
           
             
               
                 
                   
                     cov 
                     ⁡ 
                     
                       ( 
                       
                         
                           
                             - 
                             2 
                           
                           ⁢ 
                           
                             lnp 
                             i 
                           
                         
                         , 
                         
                           
                             - 
                             2 
                           
                           ⁢ 
                           l 
                           ⁢ 
                           n 
                           ⁢ 
                           
                             p 
                             j 
                           
                         
                       
                       ) 
                     
                   
                   = 
                   
                     { 
                     
                       
                         
                           
                             
                               
                                 ρ 
                                 
                                   i 
                                   ⁢ 
                                   j 
                                 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     
                                       3 
                                       . 
                                       2 
                                     
                                     ⁢ 
                                     5 
                                   
                                   + 
                                   
                                     
                                       0 
                                       . 
                                       7 
                                     
                                     ⁢ 
                                     5 
                                     ⁢ 
                                     
                                       ρ 
                                       
                                         i 
                                         ⁢ 
                                         j 
                                       
                                     
                                   
                                 
                                 ) 
                               
                             
                             , 
                           
                         
                         
                           
                             0 
                             ≤ 
                             
                               ρ 
                               
                                 i 
                                 ⁢ 
                                 j 
                               
                             
                             ≤ 
                             1 
                           
                         
                       
                       
                         
                           
                             
                               
                                 ρ 
                                 
                                   i 
                                   ⁢ 
                                   j 
                                 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     
                                       3 
                                       . 
                                       2 
                                     
                                     ⁢ 
                                     7 
                                   
                                   + 
                                   
                                     
                                       0 
                                       . 
                                       7 
                                     
                                     ⁢ 
                                     1 
                                     ⁢ 
                                     
                                       ρ 
                                       
                                         i 
                                         ⁢ 
                                         j 
                                       
                                     
                                   
                                 
                                 ) 
                               
                             
                             , 
                           
                         
                         
                           
                             
                               - 
                               0.5 
                             
                             ≤ 
                             
                               ρ 
                               
                                 i 
                                 ⁢ 
                                 j 
                               
                             
                             &lt; 
                             0 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   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 χ 2  distribution, where T=T 0 /c with 2N/c degrees of freedom and c=σ 2 /4N. 
     Here, we test if the expression levels of the subset of genes corresponding to the cellular and molecular entities used in our model align with a given model-predicted discrete state profile. Our null hypothesis is that the experimental measures do not align with model predictions of greater than control (normal), lower than the control or in alignment with the control levels. First, the p-values for individual variables, p i , were calculated using two-sample t-tests comparing the expression levels in ill subjects with those measured in healthy controls. Where the model predicts marker expression to be high (+1), a ‘right-handed’ one-tailed test was used to confirm that the measured expression levels are significantly greater than the reference control (0). Conversely, a ‘left-handed’ test was used when the model predicted a low (−1), to confirm that the measured expression levels are, on average, significantly lower than the reference control (0). For the case where the model predicted normal expression levels for a variable (0), a two-tailed t-test was used. However, the p-value from the two-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 took 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. The unadulterated correlation values ρ ij  between two variables i and j were calculated in the control groups as the pairwise Pearson&#39;s linear correlation coefficient between variables. Aggregate p-values between the predicted model states were determined via Brown&#39;s method above, where values between the predicted states that were found to disagree with the p-values for individual variables, p i , were taken as 1. Conversely, when the values between the predicted states were found to agree, p i  was assigned a standard minimum value for significance of 0.05 to avoid numerical instability in the calculation of the combined χ 2  statistic T 0 . Hold-out cross-validation was performed by randomly splitting the healthy control dataset evenly in two and performing the above method 100 times. This resulted in a maximum of 94% of the runs, finding that the test group aligned with the healthy SS0 stable state with an even distribution of the remaining cases being assigned to the other alternate stable states. 
     A comparison of the measured states with the model-predicted stable states is best visualized by projecting the multi-dimensional co-expression profiles into a two-dimensional space using multidimensional scaling. Here, the dissimilarity matrix defined by the p-value is scaled such that the 2D Euclidean distances between points approximate of the corresponding dissimilarities. This is performed using the function mdscale in MATLAB to minimize Kruskal&#39;s stress criterion normalized by the sum of squares of the dissimilarities. To maintain the relative structure of the healthy control baseline-defined stable state positions between the conditions, the weight of the added illness condition distances was set at ten percent of the distances defining the stable state positions. This value was chosen to keep the RMSD of the stable states between conditions equivalent to two decimal places. The resulting 2D plot illustrates the statistical significance of separation between the measured and predicted co-expression patterns. To compare the plots between the conditions, all the resulting 2D plots were aligned with the human baseline condition via an affine geometric transformation based on the model-predicted stable state positions via the MATLAB function fitgeotrans. 
     To provide a measure of agreement between conditions, the root mean square deviation (RMSD) between the aligned 2D plots was calculated as 
     
       
         
           
             
               
                 
                   RMSD 
                   = 
                   
                     
                       
                         
                           1 
                           N 
                         
                         ⁢ 
                         
                           
                             ∑ 
                             
                               i 
                               = 
                               1 
                             
                             N 
                           
                           ⁢ 
                           
                             
                               ( 
                               
                                 
                                   x 
                                   i 
                                   B 
                                 
                                 - 
                                 
                                   x 
                                   i 
                                   A 
                                 
                               
                               ) 
                             
                             2 
                           
                         
                       
                       + 
                       
                         
                           ( 
                           
                             
                               y 
                               i 
                               B 
                             
                             - 
                             
                               y 
                               i 
                               A 
                             
                           
                           ) 
                         
                         2 
                       
                     
                     2 
                   
                 
               
               
                 
                   ( 
                   6 
                   ) 
                 
               
             
           
         
       
     
     where x i   j  and y i   j  are the horizontal and vertical positions of the i th  point of n total points in plot j. 
     9. Simulating Intervention Courses 
     To identify a robust sequence of interventions capable of moving the system from a pathological mode of regulation to that of normal health, we evolved solutions combining a specific choice of treatment targets as well as the sequence, spacing, and type of external perturbation. For each of these candidate treatment courses, simulations were conducted to evaluate the occurrence of normal homeostasis. Specifically, each clinical intervention was represented as a treatment vector with n variables. Interventions applied to the system state at some point in time, t, are represented by the vector {right arrow over (T)}(t), such that 
         {right arrow over (T)} ( t )=( T   1 ( t ), T   2   , . . . ,T   N ( t ))  (7)
 
     where T i  is a ternary value describing the effect of the clinical treatment on the i th  element of the system: −1 (inhibiting), 0 (untreated), and 1 (supplementing). At those times where an intervention is being applied, the image vector {right arrow over (x)}(t+1) describing the preferred state towards which the system should evolve is now defined as 
         {right arrow over (x)} ( t+ 1)= {right arrow over (x)} ( t )+ {right arrow over (T)} ( t )  (8)
 
     as opposed to the unperturbed logic described in Equation (2). Due to the ternary nature of this system no value can extend beyond the range of −1 to 1; hence, values beyond this range were rounded accordingly (i.e., if x i (t)=1 and T i (t)=1, then x i (t+1)=x i (t)+T(t)=2 is rounded to 1). At times t, when there is no treatment, applied (i.e., all T(t)=0), the state transition continues according to the logic in Equation (2). 
     10. Genetic Algorithm for Optimizing Treatment Course 
     A first mapping of the illness basin of attraction consisted of a series of simulations where we first allowed for a one-time simultaneous perturbation of two or more variables at the outset only without subsequent intervention. We then expanded on this by performing a global search to find an optimal series of single-target interventions separated in time, which reliably led to health. A Genetic Algorithm (GA)-based search [66] was used to optimize this treatment course, as its form naturally accommodates the discrete definition of each system state. A treatment course vector {right arrow over (C)} with M interventions is therefore defined as 
         {right arrow over (C)} =({right arrow over ( T   1 )}( t   1start   ,t   1stop ),{right arrow over ( T   2 )}( t   2start   ,t   2stop ), . . . ,{right arrow over ( T   M )}( t   Mstart   ,t   Mstop ))  (9)
 
     where {right arrow over (T)} j (t istart , t istop ) is the j th  intervention treatment vector starting at time point t jstart  and ending at t jstop . Due to the asynchronous nature of the model each treatment vector {right arrow over (T)} only contains a single target intervention T i  that affects the i th  element variable in the system at any given time step. 
     The GA starts by generating a population of 1000 candidate treatment courses each composed of a specific number of randomly selected interventions applied at random time points. The response of the system to each treatment course in this initial generation of candidates is then simulated for 1000 timesteps. Over the course of these time steps the state of the system evolves according to Equation (2), except at those times when interventions are applied. At these intervention events the state transition follows Equation (5). These 1000 iterations provide a distribution of paths that are then ranked according to a fitness function based on the number of times a treatment successfully reaches the healthy stable state (% HHM). After all treatments in the generation have been ranked, the top 10% are retained without change for the next generation. Thus, treatments generated by the GA that have a higher probability of leading to the HEIM are re-executed and re-evaluated over many generations. The remainder of the next generation of candidate solutions is created by choosing random pairs from the total set of treatment courses (including the top 10 th  percentile) and combining (cross-over recombination) them. Combination ⊗, of two treatment courses {right arrow over (C)} 1 , and {right arrow over (C)} 2  is performed at a single random splice timepoint s to create a new treatment course {right arrow over (C)}′ 1  for the next generation, such that all treatments in {right arrow over (C)} 1  preceding s, and all treatments in {right arrow over (C)} 2  after s are included in {right arrow over (C)}′ 1 . 
     The response to the new treatment courses, including members of the previous top 10th percentile, are then simulated once again and ranked as a new generation. This process is continued iteratively for 1000 generations. The final treatment course with the highest % HEM was taken as the best treatment solution for a given run. The overall best treatment course was chosen from the 100 repeated GA runs. 
     5. Conclusion 
     Our simulations predict treatment avenues for GWI in totality, as well as specific treatments for GWI with probable PTSD. However, it must be noted that the predictions presented here are based on an idealized model of immune interactions across the BBB. The alignment of experimental data with model predictions shows that GWI aligns with a state characterized by BBB dysfunction as predicted by the model; however, additional molecular and biochemical evidence (e.g., small molecule leakage, albumin leakage, MMP activation, and monocyte adhesion to the microvascular endothelium) could be beneficial to make a more definitive claim concerning BBB function in GWI and GWI mouse models. As such, future research in this area may be useful. Overall, the advantage of this modeling methodology may also be considered one of its limitations in that it does not account for detailed kinetics. The magnitude and transition time of the interactions between elements of the extended neuroendocrine-immune system are not predicted. Rather only the logical feasibility of these interventions is assessed. While refinement of the simulation parameters based on detailed measurement of the model elements can serve to improve the accuracy of the simulation timing and magnitude, the underlying logic is not affected. Regardless, even with these refinements, the safety and efficacy of any predicted treatment strategy should be determined clinically. 
     Example Computer System 
       FIG. 13  illustrates a block diagram of an exemplary system for managing and communicating the process steps above. The system shown in  FIG. 13  is only one example of a suitable system and is not intended to limit the scope of use or functionality of embodiments of the present invention described above. The system is operational with numerous other general purpose or special purpose computing system environments or configurations. Examples of well-known computing systems, environments, and/or configurations that may be suitable for use with the information processing system include, but are not limited to, personal computer systems, server computer systems, thin clients, thick clients, hand-held or laptop devices, multiprocessor systems, microprocessor-based systems, set-top boxes, programmable consumer electronics, network PCs, minicomputer systems, mainframe computer systems, clusters, and distributed cloud computing environments that include any of the above systems or devices, and the like. 
     The system may be described in the general context of computer-executable instructions being executed by a computer system. The system may be practiced in various computing environments such as conventional and distributed computing environments where tasks are performed by remote processing devices that are linked through a communications network. In a distributed computing environment, program modules may be located in both local and remote computer system storage media, including memory storage devices. 
     Referring again to  FIG. 13 , the system includes server  700 . In some embodiments, server  700  can be embodied as a general-purpose computing device. The components of server  700  can include, but are not limited to, one or more processor devices or processing units  704 , a system memory  706 , and a bus  708  that couples various system components, including the system memory  706  to the processor  704 . 
     The bus  708  represents one or more of any of several types of bus structures, including a memory bus or memory controller, a peripheral bus, an accelerated graphics port, and a processor or local bus using any of a variety of bus architectures. By way of example, and not limitation, such architectures include Industry Standard Architecture (ISA) bus, Micro Channel Architecture (MCA) bus, Enhanced ISA (EISA) bus, Video Electronics Standards Association (VESA) local bus, and Peripheral Component Interconnects (PCI) bus. 
     The system memory  706  can also include computer system readable media in the form of volatile memory, such as random access memory (RAM)  710  and/or cache memory  712 . The server  700  can further include other removable/non-removable, volatile/non-volatile computer system storage media. By way of example only, a storage system  714  can be provided for reading from and writing to a non-removable or removable, non-volatile media such as one or more solid state disks and/or magnetic media (typically called a “hard drive”). A magnetic disk drive for reading from and writing to a removable, non-volatile magnetic disk (e.g., a “floppy disk”), and an optical disk drive for reading from or writing to a removable, non-volatile optical disk such as a CD-ROM, DVD-ROM or other optical media can be provided. In such instances, each can be connected to the bus  708  by one or more data media interfaces. The memory  706  can include at least one program product embodying a set of program modules  718  that are configured to carry out one or more features and/or functions of the present invention e.g., described with reference to  FIGS. 1-12 . Referring again to  FIG. 13 , program/utility  716 , having a set of program modules  718 , may be stored in memory  706  by way of example, and not limitation, as well as an operating system, one or more application programs, other program modules, and program data. Generally, program modules may include routines, programs, objects, components, logic, data structures, and so on that perform particular tasks or implement particular abstract data types. Each of the operating system, one or more application programs, other program modules, and program data or some combination thereof, may include an implementation of a networking environment. In some embodiments, program modules  1018  are configured to carry out one or more functions and/or methodologies of embodiments of the present invention. 
     The server  700  can also communicate with one or more external devices  1020  that enable interaction with the server  700 ; and/or any devices (e.g., network card, modem, etc.) that enable communication with one or more other computing devices. A few (non-limiting) examples of such devices include: a keyboard, a pointing device, a display  722  presenting system performance tuning dashboard  500 , etc.; one or more devices that enable a user to interact with the server  700 ; and/or any devices (e.g., network card, modem, etc.) that enable the server  700  to communicate with one or more other computing devices. Such communication can occur via I/O interfaces  724 . In some embodiments, the server  700  can communicate with one or more networks such as a local area network (LAN), a general wide area network (WAN), and/or a public network (e.g., the Internet) via network adapter  726 , enabling the system  700  to access parameters (not shown). As depicted, the network adapter  726  communicates with the other components of the server  700  via the bus  708 . Other hardware and/or software components can also be used in conjunction with the server  700 . Examples include but are not limited to: microcode, device drivers, redundant processing units, external disk drive arrays, RAID systems, tape drives, and data archival storage systems. 
     As will be appreciated by one skilled in the art, aspects of the present invention may be embodied as a system, method, or computer program product  790  at any possible technical detail level of integration. The computer program product  790  may include a computer readable storage medium (or media) having computer readable program instructions thereon for causing a processor to carry out aspects of the present invention. 
     The computer readable storage medium can be a tangible device that can retain and store instructions for use by an instruction execution device. The computer readable storage medium may be, for example, but is not limited to, an electronic storage device, a magnetic storage device, an optical storage device, an electromagnetic storage device, a semiconductor storage device, or any suitable combination of the foregoing. A non-exhaustive list of more specific examples of the computer readable storage medium includes the following: a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), a static random access memory (SRAM), a portable compact disc read-only memory (CD-ROM), a digital versatile disk (DVD), a memory stick, a floppy disk, a mechanically encoded device such as punch-cards or raised structures in a groove having instructions recorded thereon, and any suitable combination of the foregoing. A computer readable storage medium, as used herein, is not to be construed as being transitory signals per se, such as radio waves or other freely propagating electromagnetic waves, electromagnetic waves propagating through a waveguide or other transmission media (e.g., light pulses passing through a fiber-optic cable), or electrical signals transmitted through a wire. 
     Accordingly, aspects of the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment (including firmware, resident software, micro-code, etc.) or an embodiment combining software and hardware aspects that may all generally be referred to herein as a “circuit,” “module” or “system.” Furthermore, aspects of the present invention may take the form of a computer program product embodied in one or more computer readable medium(s) having computer readable program code embodied thereon. 
     Any combination of one or more computer readable medium(s) may be utilized. The computer readable medium may be a computer readable signal medium or a computer readable storage medium. A computer readable storage medium may be, for example, although not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of the computer readable storage medium would include the following: an electrical connection having one or more wires, a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing. In the context of this document, a computer readable storage medium may be any tangible medium that can contain or store a program for use by or in connection with an instruction execution system, apparatus, or device. 
     A computer readable signal medium may include a propagated data signal with computer readable program code embodied therein, for example, in baseband or as part of a carrier wave. Such a propagated signal may take any of a variety of forms, including, although not limited to, electro-magnetic, optical, or any suitable combination thereof. A computer readable signal medium may be any computer readable medium that is not a computer readable storage medium and that can communicate, propagate, or transport a program for use by or in connection with an instruction execution system, apparatus, or device. 
     Computer readable program instructions described herein can be downloaded to respective computing/processing devices from a computer readable storage medium or to an external computer or external storage device via a network, for example, the Internet, a local area network, a wide area network and/or a wireless network. The network may comprise copper transmission cables, optical transmission fibers, wireless transmission, routers, firewalls, switches, gateway computers and/or edge servers. A network adapter card or network interface in each computing/processing device receives computer readable program instructions from the network and forwards the computer readable program instructions for storage in a computer readable storage medium within the respective computing/processing device. 
     Computer readable program instructions for carrying out operations of the present invention may be assembler instructions, instruction-set-architecture (ISA) instructions, machine instructions, machine dependent instructions, microcode, firmware instructions, state-setting data, configuration data for integrated circuitry, or either source code or object code written in any combination of one or more programming languages, including an object oriented programming language such as Smalltalk, C++, or the like, and procedural programming languages, such as the “C” programming language or similar programming languages. The computer readable program instructions may execute entirely on the user&#39;s computer, partly on the user&#39;s computer, as a stand-alone software package, partly on the user&#39;s computer and partly on a remote computer or entirely on the remote computer or server. In the latter scenario, the remote computer may be connected to the user&#39;s computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection may be made to an external computer (for example, through the Internet using an Internet Service Provider). In some embodiments, electronic circuitry including, for example, programmable logic circuitry, field-programmable gate arrays (FPGA), or programmable logic arrays (PLA) may execute the computer readable program instructions by utilizing state information of the computer readable program instructions to personalize the electronic circuitry, in order to perform aspects of the present invention. 
     Aspects of the present invention have been discussed above with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems) and computer program products according to various embodiments of the invention. It will be understood that each block of the flowchart illustrations and/or block diagrams, and combinations of blocks in the flowchart illustrations and/or block diagrams, can be implemented by computer program instructions. 
     These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks. These computer program instructions may also be stored in a non-transitory computer readable storage medium that can direct a computer, other programmable data processing apparatus, or other devices to function in a particular manner, such that the instructions stored in the computer readable medium produce an article of manufacture including instructions which implement the function/act specified in the flowchart and/or block diagram block or blocks. 
     The computer program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other devices to cause a series of operational steps to be performed on the computer, other programmable apparatus or other devices to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide processes for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks. 
     The flowchart and block diagrams in the Figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of the present invention. In this regard, each block in the flowchart or block diagrams may represent a module, segment, or portion of instructions, which comprises one or more executable instructions for implementing the specified logical function(s). In some alternative implementations, the functions noted in the blocks may occur out of the order noted in the Figures. For example, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams and/or flowchart illustration, and combinations of blocks in the block diagrams and/or flowchart illustration, can be implemented by special purpose hardware-based systems that perform the specified functions or acts or carry out combinations of special purpose hardware and computer instructions. 
     The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. As used herein, the singular forms “a”, “an” and “the” are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms “comprises” and/or “comprising,” when used in this specification, specify the presence of stated features, integers, steps, operations, elements, and/or components, although do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and/or groups thereof. 
     The description of the present application has been presented for purposes of illustration and description, but is not intended to be exhaustive or limited to the invention in the form disclosed. Many modifications and variations will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the invention. The embodiments were chosen and described in order to best explain the principles of the invention and the practical application, and to enable others of ordinary skill in the art to understand various embodiments of the present invention, with various modifications as are suited to the particular use contemplated. 
     LIST OF CITED AND INCORPORATED REFERENCES 
     
         
         1. Binns, J.; Bloom, F.; Bunker, J.; Crawford, F.; Golomb, B.; Graves, J.; Klimas, N.; O&#39;&#39;Callaghan, J.; Ondra, S.; Philbert, M.; et al. Research advisory committee on Gulf War veterans&#39; illnesses.  Gulf War illness and the health of gulf war veterans: Research update and recommendations,  2009-2013 . Washington, D.C.: U.S. Government Printing Office  2014. 
         2. Golier, J. A.; Schmeidler, J.; Legge, J.; Yehuda, R. Twenty-four hour plasma cortisol and adrenocorticotropic hormone in Gulf War veterans: relationships to posttraumatic stress disorder and health symptoms.  Biological Psychiatry  2007, 62, 1175-1178. 
         3. Golier, J. A.; Caramanica, K.; Michaelides, A. C.; Makotkine, I.; Schmeidler, J.; Harvey, P. D.; Yehuda, R. A randomized, double-blind, placebo-controlled, crossover trial of mifepristone in Gulf War veterans with chronic multisymptom illness.  Psychoneuroendocrinology  2016, 64, 22-30. 
         4. Golier, J. A.; Schmeidler, J.; Yehuda, R. Pituitary response to metyrapone in Gulf War veterans: relationship to deployment, PTSD and unexplained health symptoms.  Psychoneuroendocrinology  2009, 34, 1338-1345. 
         5. Golier, J. A.; Schmeidler, J.; Legge, J.; Yehuda, R. Enhanced cortisol suppression to dexamethasone associated with Gulf War deployment.  Psychoneuroendocrinology  2006, 31, 1181-1189. 
         6. Broderick, G.; Ben-Hamo, R.; Vashishtha, S.; Efroni, S.; Nathanson, L.; Barnes, Z.; Fletcher, M. A.; Klimas, N. Altered immune pathway activity under exercise challenge in Gulf War Illness: an exploratory analysis.  Brain, behavior, and immunity  2013, 28, 159-169. 
         7. Janulewicz, P.; Krengel, M.; Quinn, E.; Heeren, T.; Toomey, R.; Killiany, R.; Zundel, C.; Ajama, J.; O&#39;Callaghan, J.; Steele, L. The multiple hit hypothesis for gulf war illness: Self-reported chemical/biological weapons exposure and mild traumatic brain injury.  Brain sciences  2018, 8, 198. 
         8. Amourette, C.; Lamproglou, I.; Barbier, L.; Fauquette, W.; Zoppe, A.; Viret, R.; Diserbo, M. Gulf War illness: Effects of repeated stress and pyridostigmine treatment on blood-brain barrier permeability and cholinesterase activity in rat brain.  Behavioural brain research  2009, 203, 207-214. 
         9. Lamproglou, I.; Barbier, L.; Diserbo, M.; Fauvelle, F.; Fauquette, W.; Amourette, C. Repeated stress in combination with pyridostigmine: Part I: Long-term behavioural consequences.  Behavioural brain research  2009, 197, 301-310. 
         10. Barbier, L.; Diserbo, M.; Lamproglou, I.; Amourette, C.; Peinnequin, A.; Fauquette, W. Repeated stress in combination with pyridostigmine: part II: changes in cerebral gene expression.  Behavioural brain research  2009, 197, 292-300. 
         11. Georgopoulos, A. P.; James, L. M.; Mahan, M. Y.; Joseph, J.; Georgopoulos, A.; Engdahl, B. E. Reduced human leukocyte antigen (HLA) protection in Gulf War Illness (GWI).  EBioMedicine  2016, 3, 79-85. 
         12. Radhakrishnan, K.; Hauser, E. R.; Polimanti, R.; Helmer, D. A.; Provenzale, D.; McNeil, R. B.; Maffucci, A.; Quaden, R.; Zhao, H.; Whitbourne, S. B. Genomics of Gulf War Illness in US Veterans Who Served during the 1990-1991 Persian Gulf War: Methods and Rationale for Veterans Affairs Cooperative Study #2006 . Brain sciences  2021, 11, 845. 
         13. Alshelh, Z.; Albrecht, D. S.; Bergan, C.; Akeju, O.; Clauw, D. J.; Conboy, L.; Edwards, R. R.; Kim, M.; Lee, Y. C.; Protsenko, E. In-vivo imaging of neuroinflammation in veterans with Gulf War illness.  Brain, behavior, and immunity  2020, 87, 498-507. 
         14. Engdahl, B. E.; James, L. M.; Miller, R. D.; Leuthold, A. C.; Lewis, S. M.; Carpenter, A. F.; Georgopoulos, A. P. A magnetoencephalographic (MEG) study of Gulf War Illness (GWI).  EBioMedicine  2016, 12, 127-132. 
         15. Rayhan, R. U.; Stevens, B. W.; Raksit, M. P.; Ripple, J. A.; Timbol, C. R.; Adewuyi, O.; VanMeter, J. W.; Baraniuk, J. N. Exercise challenge in Gulf War Illness reveals two subgroups with altered brain structure and function.  PLoS One  2013, 8, e63903. 
         16. Nicolson, G. L. Gulf War Illnesses: Causes and Treatments.  Armed Forces Medical Developments  2001, 2, 41-44. 
         17. Nicolson, G. L.; Nasralla, M. Y.; Franco, R.; Nicolson, N. L.; Erwin, R.; Ngwenya, R.; Berns, P. A. Diagnosis and Integrative Treatment of Intracellular Bacterial Infections in Chronic Fatigue and Fibromyalgia Syndromes, Gulf War Illness, Rheumatoid Arthritis, and Other Chronic Illnesses.  Clinical Practice of Alternative Medicine  2000, 1, 92-106. 
         18. Weiner, M. W.; Meyerhoff, D. J.; Neylan, T. C.; Hlavin, J.; Ramage, E. R.; McCoy, D.; Studholme, C.; Cardenas, V.; Marmar, C.; Truran, D. The relationship between Gulf War illness, brain N-acetylaspartate, and post-traumatic stress disorder.  Military medicine  2011, 176, 896-902. 
         19. Jeffrey, M.; Collado, F.; Kibler, J.; DeLucia, C.; Messer, S.; Klimas, N.; Craddock, T. J. Post-traumatic stress impact on health outcomes in Gulf War Illness.  BMC psychology  2021, 9, 1-9. 
         20. Jeffrey, M. G.; Krengel, M.; Kibler, J. L.; Zundel, C.; Klimas, N. G.; Sullivan, K.; Craddock, T. J. Neuropsychological findings in Gulf War illness: a review.  Frontiers in psychology  2019, 7, 2088. 
         21. Plantinga, L.; Bremner, J. D.; Miller, A. H.; Jones, D. P.; Veledar, E.; Goldberg, J.; Vaccarino, V. Association between posttraumatic stress disorder and inflammation: a twin study.  Brain, behavior, and immunity  2013, 30, 125-132. 
         22. Hoge, E.; Brandstetter, K.; Moshier, S.; Pollack, M.; Wong, K.; Simon, N. Broad spectrum of cytokine abnormalities in panic disorder and posttraumatic stress disorder.  Depression and anxiety  2009, 26, 447-455. 
         23. Guo, M.; Liu, T.; Guo, J.-C.; Jiang, X.-L.; Chen, F.; Gao, Y.-S. Study on serum cytokine levels in posttraumatic stress disorder patients.  Asian Pacific journal of tropical medicine  2012, 5, 323-325. 
         24. Baker, D. G.; Nievergelt, C. M.; O&#39;Connor, D. T. Biomarkers of PTSD: neuropeptides and immune signaling.  Neuropharmacology  2012, 62, 663-673. 
         25. Khaiboullina, S. F.; DeMeirleir, K. L.; Rawat, S.; Berk, G. S.; Gaynor-Berk, R. S.; Mijatovic, T.; Blatt, N.; Rizvanov, A. A.; Young, S. G.; Lombardi, V. C. Cytokine expression provides clues to the pathophysiology of Gulf War illness and myalgic encephalomyelitis.  Cytokine  2015, 72, 1-8. 
         26. Broderick, G.; Kreitz, A.; Fuite, J.; Fletcher, M. A.; Vernon, S. D.; Klimas, N. A pilot study of immune network remodeling under challenge in Gulf War Illness.  Brain, behavior, and immunity  2011, 25, 302-313. 
         27. Friedman, A.; Kaufer, D.; Shemer, J.; Hendler, I.; Soreq, H.; Tur-Kaspa, I. Pyridostigmine brain penetration under stress enhances neuronal excitability and induces early immediate transcriptional response.  Nature medicine  1996, 2, 1382-1385. 
         28. Esposito, P.; Chandler, N.; Kandere, K.; Basu, S.; Jacobson, S.; Connolly, R.; Tutor, D.; Theoharides, T. C. Corticotropin-releasing hormone and brain mast cells regulate blood-brain-barrier permeability induced by acute stress.  Journal of Pharmacology and Experimental Therapeutics  2002, 303, 761-1066. 
         29. Groeneweg, F. L.; Karst, H.; de Kloet, E. R.; Joëls, M. Rapid non-genomic effects of corticosteroids through the membrane-associated MR and GR and their role in the central stress.  J Endocrinol  2011, 2, 153-167. 
         30. Lupien, S. J.; McEwen, B. S.; Gunnar, M. R.; Heim, C. Effects of stress throughout the lifespan on the brain, behaviour and cognition.  Nature reviews neuroscience  2009, 7, 434-445. 
         31. Mikics, E.; Kruk, M. R.; Haller, J. Genomic and non-genomic effects of glucocorticoids on aggressive behavior in male rats.  Psychoneuroendocrinology  2004, 29, 618-635. 
         32. Rankin, J.; Walker, J. J.; Windle, R.; Lightman, S. L.; Terry, J. R. Characterizing dynamic interactions between ultradian glucocorticoid rhythmicity and acute stress using the phase response curve.  PloS one  2012, 7, e30978. 
         33. Walker, J. J.; Spiga, F.; Waite, E.; Zhao, Z.; Kershaw, Y.; Terry, J. R.; Lightman, S. L. The origin of glucocorticoid hormone oscillations.  PLoS Biol  2012, 7, e1001341. 
         34. Walker, J. J.; Terry, J. R.; Lightman, S. L. Origin of ultradian pulsatility in the hypothalamic-pituitary-adrenal axis.  Proceedings of the Royal Society B: Biological Sciences  2010, 277, 1627-1633. 
         35. Vinther, F.; Andersen, M.; Ottesen, J. T. The minimal model of the hypothalamic-pituitary-adrenal axis.  Journal of mathematical biology  2011, 63, 663-690. 
         36. Ben-Zvi, A.; Vernon, S. D.; Broderick, G. Model-based therapeutic correction of hypothalamic-pituitary-adrenal axis dysfunction.  PLoS computational biology  2009, 5, e1000273. 
         37. Gupta, S.; Aslakson, E.; Gurbaxani, B. M.; Vernon, S. D. Inclusion of the glucocorticoid receptor in a hypothalamic pituitary adrenal axis model reveals bistability.  Theoretical Biology and Medical Modelling  2007, 4, 8. 
         38. Craddock, T. J. A.; Fritsch, P.; Rice Jr, M. A.; Del Rosario, R. M.; Miller, D. B.; Fletcher, M. A.; Klimas, N. G.; Broderick, G. A role for homeostatic drive in the perpetuation of complex chronic illness: Gulf War Illness and chronic fatigue syndrome.  PLoS One  2014, 9, e84839. 
         39. Fritsch, P.; Craddock, T. J. A.; del Rosario, R. M.; Rice, M. A.; Smylie, A.; Folcik, V. A.; de Vries, G.; Fletcher, M. A.; Klimas, N. G.; Broderick, G. Succumbing to the laws of attraction: Exploring the sometimes pathogenic versatility of discrete immune logic.  Systems Biomedicine  2014, 1, 179-194. 
         40. Mendoza, L.; Xenarios, I. A method for the generation of standardized qualitative dynamical systems of regulatory networks. Theoretical Biology and Medical Modelling 2006, 3, 13. 
         41. Thomas, R.; Thieffry, D.; Kaufman, M. Dynamical behaviour of biological regulatory networks—I. Biological role of feedback loops and practical use of the concept of the loop-characteristic state.  Bulletin of mathematical biology  1995, 57, 247-276. 
         42. Thomas, R.; Kaufman, M. Multistationarity, the basis of cell differentiation and memory. II. Logical analysis of regulatory networks in terms of feedback circuits. Chaos: An  Interdisciplinary Journal of Nonlinear Science  2001, 11, 180-195. 
         43. Thomas, R. Regulatory networks seen as asynchronous automata: a logical description.  Journal of theoretical biology  1991, 153, 1-23. 
         44. Whitley, D.; Garrett, D.; Watson, J.-P. Quad search and hybrid genetic algorithms. In Proceedings of the Genetic and Evolutionary Computation Conference, 2003; pp. 1469-1480. 
         45. Day, J.; Rubin, J.; Vodovotz, Y.; Chow, C. C.; Reynolds, A.; Clermont, G. A reduced mathematical model of the acute inflammatory response II. Capturing scenarios of repeated endotoxin administration.  Journal of theoretical biology  2006, 242, 237-256. 
         46. Reynolds, A.; Rubin, J.; Clermont, G.; Day, J.; Vodovotz, Y.; Ermentrout, G. B. A reduced mathematical model of the acute inflammatory response: I. Derivation of model and analysis of anti-inflammation.  Journal of theoretical biology  2006, 242, 220-236. 
         47. McDonald, S. D.; Thompson, N. L.; Stratton, K. J.; Calhoun, P. S. Diagnostic accuracy of three scoring methods for the Davidson Trauma Scale among US military Veterans.  Journal of anxiety disorders  2014, 28, 160-168. 
         48. Storey, J. D. A direct approach to false discovery rates.  Journal of the Royal Statistical Society: Series B  ( Statistical Methodology ) 2002, 64, 479-498. 
         49. White, R. F.; Steele, L.; O&#39;Callaghan, J. P.; Sullivan, K.; Binns, J. H.; Golomb, B. A.; Bloom, F. E.; Bunker, J. A.; Crawford, F.; Graves, J. C. Recent research on Gulf War illness and other health problems in veterans of the 1991 Gulf War: Effects of toxicant exposures during deployment.  Cortex  2016, 74, 449-475. 
         50. Locker, A. R.; Michalovicz, L. T.; Kelly, K. A.; Miller, J. V.; Miller, D. B.; O&#39;Callaghan, J. P. Corticosterone primes the neuroinflammatory response to Gulf War Illness-relevant organophosphates independently of acetylcholinesterase inhibition.  Journal of Neurochemistry  2017. 
         51. O&#39;Callaghan, J. P.; Kelly, K. A.; Locker, A. R.; Miller, D. B.; Lasley, S. M. Corticosterone primes the neuroinflammatory response to DFP in mice: potential animal model of Gulf War Illness.  Journal of neurochemistry  2015, 133, 708-721. 
         52. Koo, B.-B.; Michalovicz, L. T.; Calderazzo, S.; Kelly, K. A.; Sullivan, K.; Killiany, R. J.; O&#39;Callaghan, J. P. Corticosterone potentiates DFP-induced neuroinflammation and affects high-order diffusion imaging in a rat model of Gulf War Illness.  Brain, behavior, and immunity  2018, 67, 42-46. 
         53. Craddock, T. J. A.; Michalovicz, L.; Kelly, K. A.; Rice Jr, M.; Miller, D.; Klimas, N.; Morris, M.; O&#39;Callaghan, J.; Broderick, G. A Logic Model of Neuronal-Glial Interaction Suggests Altered Homeostatic Regulation in the Perpetuation of Neuroinflammation.  Frontiers in Cellular Neuroscience  2018, 12, 336. 
         54. Miller, J. V.; LeBouf, R. F.; Kelly, K. A.; Michalovicz, L. T.; Ranpara, A.; Locker, A. R.; Miller, D. B.; O&#39;Callaghan, J. P. The neuroinflammatory phenotype in a mouse model of Gulf War Illness is unrelated to brain regional levels of acetylcholine as measured by quantitative HILIC-UPLC-MS/MS.  Toxicological Sciences  2018, 165, 302-313. 
         55. Jones, B. C.; Miller, D. B.; Lu, L.; Zhao, W.; Ashbrook, D. G.; Xu, F.; Mulligan, M. K.; Williams, R. W.; Zhuang, D.; Torres-Rojas, C. Modeling the genetic basis of individual differences in susceptibility to Gulf War Illness.  Brain sciences  2020, 7, 143. 
         56. Xu, F.; Ashbrook, D. G.; Gao, J.; Starlard-Davenport, A.; Zhao, W.; Miller, D. B.; O&#39;Callaghan, J. P.; Williams, R. W.; Jones, B. C.; Lu, L. Genome-wide transcriptome architecture in a mouse model of Gulf War Illness.  Brain, behavior, and immunity  2020, 89, 209-223. 
         57. Ashbrook, D. G.; Hing, B.; Michalovicz, L. T.; Kelly, K. A.; Miller, J. V.; de Vega, W. C.; Miller, D. B.; Broderick, G.; O&#39;Callaghan, J. P.; McGowan, P. O. Epigenetic impacts of stress priming of the neuroinflammatory response to sarin surrogate in mice: a model of Gulf War illness.  Journal of neuroinflammation  2018, 15, 1-15. 
         58. Gao, J.; Xu, F.; Starlard-Davenport, A.; Miller, D. B.; O&#39;Callaghan, J. P.; Jones, B. C.; Lu, L. Exploring the role of chemokine receptor 6 (Ccr6) in the BXD mouse model of Gulf War illness.  Frontiers in Neuroscience  2020, 14, 818. 
         59. Belgrad, J.; Dutta, D. J.; Bromley-Coolidge, S.; Kelly, K. A.; Michalovicz, L. T.; Sullivan, K. A.; O&#39;Callaghan, J. P.; Fields, R. D. Oligodendrocyte involvement in gulf war illness.  Glia  2019, 67, 2107-2124. 
         60. Lacagnina, M. J.; Li, J.; Lorca, S.; Rice, K. C.; Sullivan, K.; O&#39;Callaghan, J. P.; Grace, P. M. A role for neuroimmune signaling in a rat model of Gulf War Illness-related pain.  Brain, behavior, and immunity  2021, 91, 418-428. 
         61. Michalovicz, L. T.; Locker, A. R.; Kelly, K. A.; Miller, J. V.; Barnes, Z.; Fletcher, M. A.; Miller, D. B.; Klimas, N. G.; Morris, M.; Lasley, S. M. Corticosterone and pyridostigmine/DEET exposure attenuate peripheral cytokine expression: Supporting a dominant role for neuroinflammation in a mouse model of Gulf War Illness.  Neurotoxicology  2019, 70, 26-32. 
         62. Cheng, C.-H.; Koo, B.-B.; Calderazzo, S.; Quinn, E.; Aenlle, K.; Steele, L.; Klimas, N.; Krengel, M.; Janulewicz, P.; Toomey, R. Alterations in high-order diffusion imaging in veterans with Gulf War Illness is associated with chemical weapons exposure and mild traumatic brain injury.  Brain, behavior, and immunity  2020, 89, 281-290. 
         63. Kansagara, D.; Emily, G. G.; Ayers, C. K.; Winchell, K. A.; Freeman, M.; Press, A. M.; Paynter, R.; Kansagara, M. D.; Nugent, M. S. M. Biological Measures and Diagnostic Tools for Gulf War Illness: A Systematic Review. 2020. 
         64. Gean, E. G.; Ayers, C. K.; Winchell, K. A.; Freeman, M.; Press, A. M.; Paynter, R.; Kansagara, D.; Nugent, S. M. Biological measures and diagnostic tools for Gulf War Illness-A systematic review.  Life sciences  2021, 119360. 
         65. National Library of Medicine, U.S. TNF and Glucocorticoid Antagonist for GWI Associated Multi-symptom Disease Homeostasis Reset (E/M). Available online: https://clinicaltrials.gov/ct2/show/NCT04254627 (accessed on 5 Feb. 2020). 
         66. Craddock, T. J. A.; Del Rosario, R. R.; Rice, M.; Zysman, J. P.; Fletcher, M. A.; Klimas, N. G.; Broderick, G. Achieving Remission in Gulf War Illness: A Simulation-Based Approach to Treatment Design.  PloS one  2015, 7, e0132774. 
         67. Ishøy, T.; Andersson, A.-M.; Suadicani, P.; Guldager, B.; Appleyard, M.; Gyntelberg, F.; Skakkebaek, N. Major reproductive health characteristics in male Gulf War Veterans. The Danish Gulf War Study.  Danish Medical Bulletin  2001, 48, 29-32. 
         68. Ishoy, T.; Suadicani, P.; Andersson, A.-M.; Guldager, B.; Appleyard, M.; Skakkebæk, N.; Gyntelberg, F. Prevalence of male sexual problems in the Danish Gulf War Study.  Scandinavian Journal of Sexology  2001. 
         69. Ishøy, T.; Suadicani, P.; Guldager, B.; Appleyard, M.; Hein, H. O.; Gyntelberg, F. State of health after deployment in the Persian Gulf. The Danish Gulf war study.  Danish medical bulletin  1999, 46, 416-419. 
         70. McDiarmid, M. A.; Squibb, K.; Engelhardt, S.; Oliver, M.; Gucer, P.; Wilson, P. D.; Kane, R.; Kabat, M.; Kaup, B.; Anderson, L. Surveillance of depleted uranium exposed Gulf War veterans: health effects observed in an enlarged “friendly fire” cohort.  Journal of occupational and environmental medicine  2001, 43, 991-1000. 
         71. McDiarmid, M.; Engelhardt, S.; Dorsey, C.; Oliver, M.; Gucer, P.; Gaitens, J.; Kane, R.; Cernich, A.; Kaup, B.; Hoover, D. Longitudinal health surveillance in a cohort of Gulf War veterans 18 years after first exposure to depleted uranium.  Journal of Toxicology and Environmental Health, Part A  2011, 74, 678-691. 
         72. McDiarmid, M. A.; Engelhardt, S. M.; Oliver, M.; Gucer, P.; Wilson, P. D.; Kane, R.; Cernich, A.; Kaup, B.; Anderson, L.; Hoover, D. Health surveillance of Gulf War I veterans exposed to depleted uranium: updating the cohort.  Health physics  2007, 93, 60-73. 
         73. McDiarmid, M. A.; Engelhardt, S.; Oliver, M.; Gucer, P.; Wilson, P. D.; Kane, R.; Kabat, M.; Kaup, B.; Anderson, L.; Hoover, D. Health effects of depleted uranium on exposed Gulf War veterans: a 7-year follow-up.  Journal of Toxicology and Environmental Health, Part A  2004, 67, 277-296. 
         74. McDiarmid, M. A.; Keogh, J. P.; Hooper, F. J.; McPhaul, K.; Squibb, K.; Kane, R.; DiPino, R.; Kabat, M.; Kaup, B.; Anderson, L. Health effects of depleted uranium on exposed Gulf War veterans.  Environmental research  2000, 82, 168-180. 
         75. McDiarmid, M. A.; Gaitens, J. M.; Hines, S.; Breyer, R.; Wong-You-Cheong, J. J.; Engelhardt, S. M.; Oliver, M.; Gucer, P.; Kane, R.; Cernich, A. The Gulf War depleted uranium cohort at 20 years: bioassay results and novel approaches to fragment surveillance.  Health physics  2013, 74, 347-361. 
         76. McDiarmid, M. A.; Engelhardt, S. M.; Oliver, M.; Gucer, P.; Wilson, P. D.; Kane, R.; Kabat, M.; Kaup, B.; Anderson, L.; Hoover, D. Biological monitoring and surveillance results of Gulf War I veterans exposed to depleted uranium.  International archives of occupational and environmental health  2006, 79, 11-21. 
         77. Joshi, S. C.; Mathur, R.; Gulati, N. Testicular toxicity of chlorpyrifos (an organophosphate pesticide) in albino rat.  Toxicology and Industrial Health  2007, 23, 439-444. 
         78. Peiris, D. C.; Dhanushka, T. Low doses of chlorpyrifos interfere with spermatogenesis of rats through reduction of sex hormones.  Environmental Science and Pollution Research  2017, 24, 20859-20867. 
         79. Alaa-Eldin, E. A.; El-Shafei, D. A.; Abouhashem, N. S. Individual and combined effect of chlorpyrifos and cypermethrin on reproductive system of adult male albino rats.  Environmental Science and Pollution Research  2017, 24, 1532-1543. 
         80. Kudavidanage, E. P.; Dissanayake, D.; Keerthirathna, W.; Nishshanke, N.; Peiris, L. D. C. Commercial formulation of chlorpyrifos alters neurological behaviors and fertility.  Biology  2020, 9, 49. 
         81. Li, J.; Ren, F.; Li, Y.; Luo, J.; Pang, G. Chlorpyrifos induces metabolic disruption by altering levels of reproductive hormones.  Journal of agricultural and food chemistry  2019, 67, 7553-10562. 
         82. White, H. D.; Robinson, T. D. A novel use for testosterone to treat central sensitization of chronic pain in fibromyalgia patients.  International Immunopharmacology  2015, 27, 244-248. 
         83. White, H. D.; Brown, L. A.; Gyurik, R. J.; Manganiello, P. D.; Robinson, T. D.; Hallock, L. S.; Lewis, L. D.; Yeo, K.-T. J. Treatment of pain in fibromyalgia patients with testosterone gel: pharmacokinetics and clinical response.  International immunopharmacology  2015, 27, 249-256. 
         84. Rice Jr, M. A.; Craddock, T. J. A.; Folcik, V. A.; del Rosario, R. M.; Barnes, Z. M.; Klimas, N. G.; Fletcher, M. A.; Zysman, J.; Broderick, G. Gulf War Illness: Is there lasting damage to the endocrine-immune circuitry?  Systems Biomedicine  2014, 2, 80-89. 
         85. Maes, M.; Lin, A.-h.; Delmeire, L.; Van Gastel, A.; Kenis, G.; De Jongh, R.; Bosmans, E. Elevated serum interleukin-6 (IL-6) and IL-6 receptor concentrations in posttraumatic stress disorder following accidental man-made traumatic events.  Biological psychiatry  1999, 45, 833-839. 
         86. Lima, B. B.; Hammadah, M.; Wilmot, K.; Pearce, B. D.; Shah, A.; Levantsevych, O.; Kaseer, B.; Obideen, M.; Gafeer, M. M.; Kim, J. H. Posttraumatic stress disorder is associated with enhanced interleukin-6 response to mental stress in subjects with a recent myocardial infarction.  Brain, behavior, and immunity  2019, 75, 26-33. 
         87. Passos, I. C.; Vasconcelos-Moreno, M. P.; Costa, L. G.; Kunz, M.; Brietzke, E.; Quevedo, J.; Salum, G.; Magalhães, P. V.; Kapczinski, F.; Kauer-Sant&#39;Anna, M. Inflammatory markers in post-traumatic stress disorder: a systematic review, meta-analysis, and meta-regression.  The Lancet Psychiatry  2015, 2, 702-1012. 
         88. Cheng, Y.; Jope, R. S.; Beurel, E. A pre-conditioning stress accelerates increases in mouse plasma inflammatory cytokines induced by stress.  BMC neuroscience  2015, 16, 1-8. 
         89. Tao, W.; Dong, Y.; Su, Q.; Wang, H.; Chen, Y.; Xue, W.; Chen, C.; Xia, B.; Duan, J.; Chen, G. Liquiritigenin reverses depression-like behavior in unpredictable chronic mild stress-induced mice by regulating PI3K/Akt/mTOR mediated BDNF/TrkB pathway.  Behavioural Brain Research  2016, 308, 177-186. 
         90. Tagliari, B.; Tagliari, A. P.; Schmitz, F.; Da Cunha, A. A.; Dalmaz, C.; Wyse, A. T. Chronic variable stress alters inflammatory and cholinergic parameters in hippocampus of rats.  Neurochemical research  2011, 36, 487-493. 
         91. Reader, B. F.; Jarrett, B. L.; McKim, D. B.; Wohleb, E. S.; Godbout, J. P.; Sheridan, J. F. Peripheral and central effects of repeated social defeat stress: monocyte trafficking, microglial activation, and anxiety.  Neuroscience  2015, 289, 429-442. 
         92. Gamble-George, J. C.; Baldi, R.; Halladay, L.; Kocharian, A.; Hartley, N.; Silva, C. G.; Roberts, H.; Haymer, A.; Marnett, L. J.; Holmes, A. Cyclooxygenase-2 inhibition reduces stress-induced affective pathology.  Elife  2016, 5, e14137. 
         93. Kao, C.-Y.; He, Z.; Zannas, A. S.; Hahn, O.; Kühne, C.; Reichel, J. M.; Binder, E. B.; Wotjak, C. T.; Khaitovich, P.; Turck, C. W. Fluoxetine treatment prevents the inflammatory response in a mouse model of posttraumatic stress disorder.  Journal of psychiatric research  2016, 76, 74-83. 
         94. Ebenezer, P. J.; Wilson, C. B.; Wilson, L. D.; Nair, A. R. The anti-inflammatory effects of blueberries in an animal model of post-traumatic stress disorder (PTSD).  PloS one  2016, 11, e0160923. 
         95. Young, M. B.; Howell, L. L.; Hopkins, L.; Moshfegh, C.; Yu, Z.; Clubb, L.; Seidenberg, J.; Park, J.; Swiercz, A. P.; Marvar, P. J. A peripheral immune response to remembering trauma contributes to the maintenance of fear memory in mice.  Psychoneuroendocrinology  2018, 94, 143-151. 
         96. Fukuda, K.; Nisenbaum, R.; Stewart, G.; Thompson, W. W.; Robin, L.; Washko, R. M.; Noah, D. L.; Barrett, D. H.; Randall, B.; Herwaldt, B. L. Chronic multisymptom illness affecting Air Force veterans of the Gulf War.  Jama  1998, 280, 981-988. 
         97. Reeves, W. C.; Lloyd, A.; Vernon, S. D.; Klimas, N.; Jason, L. A.; Bleijenberg, G.; Evengard, B.; White, P. D.; Nisenbaum, R.; Unger, E. R. Identification of ambiguities in the 1994 chronic fatigue syndrome research case definition and recommendations for resolution.  BMC health services research  2003, 3, 1-9. 
         98. Collins, J. F.; Donta, S. T.; Engel Jr, C. C.; Baseman, J. B.; Dever, L. L.; Taylor, T.; Boardman, K. D.; Martin, S. E.; Wiseman, A. L.; Feussner, J. R. The antibiotic treatment trial of Gulf War Veterans&#39; Illnesses: issues, design, screening, and baseline characteristics.  Controlled clinical trials  2002, 23, 333-353. 
         99. Kelly, K. A.; Michalovicz, L. T.; Miller, J. V.; Castranova, V.; Miller, D. B.; O&#39;Callaghan, J. P. Prior exposure to corticosterone markedly enhances and prolongs the neuroinflammatory response to systemic challenge with LPS.  PloS one  2018, 13, e0190546. 
         100. Brown, M. B. 400: A method for combining non-independent, one-sided tests of significance.  Biometrics  1975, 987-992. 
       
    
     ADDITIONAL REFERENCES 
     Non-limiting representative examples demonstrating the enablement and functionality of the claimed system and process have been discussed in the following publications, which are also incorporated in their entirety by reference:
     Carrera Arias, F. J.; Aenlle, K.; Abreu, M.; Holschbach, M. A.; Michalovicz, L. T.; Kelly, K. A.; Klimas, N.; O&#39;Callaghan, J. P.; Craddock, T. J. A. Modeling Neuroimmune Interactions in Human Subjects and Animal Models to Predict Subtype-Specific Multidrug Treatments for Gulf War Illness. Int. J. Mol. Sci. 2021, 22, 8546.   

     Gene Expression Analysis to Identify Pathways And Drug Targets 
     
         
         Jeffrey, M. G., Nathanson, L., Aenlle, K., Barnes, Z. M., Baig, M., Broderick, G. &amp; Craddock, T. J. (2019). Treatment avenues in myalgic encephalomyelitis/chronic fatigue syndrome: a split-gender pharmacogenomic study of gene-expression modules. Clinical therapeutics, 41(5), 815-835. 
       
    
     Logic Modeling for Treatment Course Prediction 
     
         
         Sedghamiz, H., Morris, M., Whitley, D., Craddock, T. J., Pichichero, M., &amp; Broderick, G. (2019). Computation of robust minimal intervention sets in multi-valued biological regulatory networks. Frontiers in physiology, 10, 241. 
         Sedghamiz, H., Morris, M., Craddock, T. J., Whitley, D., &amp; Broderick, G. (2018). High-fidelity discrete modeling of the HPA axis: A study of regulatory plasticity in biology. BMC systems biology, 12(1), 76. 
         Morris, M. C., Cooney, K. E., Sedghamiz, H., Abreu, M., Collado, F., Balbin, E. G., . . . &amp; Fletcher, M. A. (2019). Leveraging prior knowledge of endocrine immune regulation in the therapeutically relevant phenotyping of women with chronic fatigue syndrome. Clinical therapeutics, 41(4), 656-674. 
         Sedghamiz, H., Morris, M., Craddock, T. J., Whitley, D., &amp; Broderick, G. (2019). Bio-ModelChecker: Using Bounded Constraint Satisfaction to Seamlessly Integrate Observed Behavior With Prior Knowledge of Biological Networks. Frontiers in Bioengineering and Biotechnology, 7, 48. 
         Lyman, C. A., Clement, M., Craddock, T. J., Fletcher, M., Klimas, N. G., &amp; Broderick, G. (2019, September). Feedback Regulation of Immune Response to Maximum Exercise in Gulf War Illness. In Proceedings of the 10th ACM International Conference on Bioinformatics, Computational Biology and Health Informatics (pp. 124-129). 
       
    
     Consensus Docking to Determine Drug Choices 
     
         
         Jaundoo, R., Bohmann, J., Gutierrez, G. E., Klimas, N., Broderick, G., &amp; Craddock, T. J. (2018). Using a consensus docking approach to predict adverse drug reactions in combination drug therapies for gulf war illness. International journal of molecular sciences, 19(11), 3355. 
         Jaundoo, R., Bohmann, J., Gutierrez, G. E., Klimas, N., Broderick, G., &amp; Craddock, T. J. (2020). Towards a Treatment for Gulf War Illness: A Consensus Docking Approach. Military Medicine, 185 (Supplement_1), 554-561. 
       
    
     Druggable Pathways Database 
     
         
         Jaundoo, R., &amp; Craddock, T. J. (2020). DRUGPATH: The Drug Gene Pathway Meta-Database. International Journal of Molecular Sciences, 21(9), 3171.