Computational Analysis of Therapeutic Neuroadaptation to Chronic Antidepressant in a Model of the Monoaminergic Neurotransmitter and Stress Hormone Systems

The clinical practice of selective serotonin reuptake inhibitor (SSRI) augmentation relies heavily on trial-and-error. Unfortunately, the drug combinations prescribed today fail to provide relief for many depressed patients. In order to identify potentially more effective treatments, we developed a computational model of the monoaminergic neurotransmitter and stress-steroid systems that neuroadapts to chronic administration of combinations of antidepressant drugs and hormones by adjusting the strengths of its transmitter-system components (TSCs). We used the model to screen 60 chronically administered drug/hormone pairs and triples, and identified as potentially therapeutic those combinations that raised the monoamines (serotonin, norepinephrine, and dopamine) but lowered cortisol following neuroadaptation in the model. We also evaluated the contributions of individual and pairs of TSCs to therapeutic neuroadaptation with chronic SSRI using sensitivity, correlation, and linear temporal-logic analyses. All three approaches revealed that therapeutic neuroadaptation to chronic SSRI is an overdetermined process that depends on multiple TSCs, providing a potential explanation for the clinical finding that no single antidepressant regimen alleviates depressive symptoms in all patients.


INTRODUCTION
Depression is a debilitating psychological disorder and a leading cause of physical disability worldwide (Cipriani et al., 2016;Friedrich, 2017). The current first-line treatment for depression is chronic administration of selective serotonin reuptake inhibitors (SSRIs), which provide complete relief from depressive symptomatology in only one-third of patients (Turner et al., 2008). In SSRI non-responders, defined as patients who experience undetectable or minimal improvement on clinical standardized rating instruments (such as the Hamilton Depression Rating Scale), clinicians employ acute or chronic alternatives to SSRI monotherapy (El-Hage et al., 2013). Acute (rapid-acting) alternatives include electroconvulsive therapy (ECT) or, more recently, ketamine, a glutamate N-methyl-D-aspartate receptor antagonist (Gideons et al., 2014). Both of these acute depression treatment approaches are promising, but current guidelines for unipolar depression treatment employ chronic antidepressant administration except in the cases of severe suicidality, and most depressed patients are still treated using chronic antidepressant administration (Davidson, 2010;Murrough et al., 2013).
Chronic alternatives include escalation of the dosage of the current SSRI, switching to a different SSRI, or augmentation with a non-SSRI antidepressant that targets the monoaminergic nervous system . Although these strategies increase the proportion of depressed patients who respond to pharmacotherapy, an effective chronic antidepressant combination unfortunately is not found in all SSRI nonresponders (Ananth, 1998;Barowsky and Schwartz, 2006).
The focus of this analysis is on chronic administration of antidepressants and other substances, either alone or in combination. There is currently no well-defined procedure for choosing which among the possible combinations of antidepressants (or combinations of antidepressants and other substances) would be most effective for a given patient. We have developed a computational model that leverages current knowledge to make predictions concerning the potential efficacy of combinations of antidepressants and the hormones mediating the hypothalamic-pituitary-adrenal (HPA) axis. The model is based on the monoaminergic neurotransmitter system and its interactions with the stress hormone system as regulated by the HPA axis.
The current approach to designing chronic antidepressant pharmocotherapy is based on the monoamine hypothesis. Mood is largely determined by the three monoaminergic neurotransmitters: serotonin (5HT), secreted from the dorsal raphe (DR); norepinephrine (NE), secreted from the locus coeruleus (LC); and dopamine (DA), secreted from the ventral tegmental area (VTA). The monoamine hypothesis was developed in the 1960s, following on clinical observations that drugs that elevate brain levels of 5HT, NE, or DA had mood-elevating effects (Schildkraut, 1965). The low efficacy-rate of SSRIs and of current SSRI augmentation strategies has led the field to investigate pharmacological targets beyond the monoamines that could improve antidepressant response efficacy, especially the interactions between the monoaminergic neurotransmitter and other neurotransmitter and hormone systems.
The most commonly observed risk factors for depression are stressful life events. For that reason, our current model is focused on the interactions between the monoaminergic neurotransmitter systems and the neuroendocrine response to stress (Kendler et al., 1999;Kendler and Gardner, 2016). Stress leads to activation of the HPA axis, resulting in elevated plasma levels of the stresssteroid, cortisol (Kendler et al., 1999). Briefly, HPA axis activation begins with corticotropin releasing factor (CRF) secretion from the paraventricular nucleus (PVN) of the hypothalamus in response to stress. Activation of pituitary gland CRF receptors by CRF results in adrenocorticotropic hormone (ACTH) release from the pituitary gland. Binding of ACTH to its receptors on the adrenal gland promotes cortisol release from the adrenal cortex (Pariante and Lightman, 2008;Morris et al., 2012). In the short-term, cortisol release in response to HPA axis activation can be beneficial for dealing with stress (Dhabhar and Mcewen, 1997;McEwen, 2004). However, chronic elevations of blood cortisol levels are related to depressive symptomology, and response to antidepressant drugs is associated with normalization of plasma cortisol levels (Johnson et al., 1992;Wong et al., 2000).
The interactions between the monoaminergic neurotransmitter system and the HPA axis are complex. For example, activation of monoaminergic receptors on the PVN, pituitary gland, and adrenal gland by monoaminergic neurotransmitters has been found to enhance HPA axis activity (Dinan, 1996;Ziegler et al., 1999;Ma and Morilak, 2005), while cortisol can alter expression of proteins involved in monoaminergic synaptic transmission, including serotonergic receptors, monoamine synthesis enzymes, and monoamine oxidase (MAO) (Hucklebridge et al., 1998;McAllister-Williams et al., 2007;Nexon et al., 2011).
Here we represent the major interactions between the monoaminergic neurotransmitter systems and the HPA axis in a computational model that we refer to as the Monoamine-Stress model (MS-model). It extends our previously published computational model of the monoaminergic neurotransmitter system but differs in its structure, training procedure, and analysis (Camacho and Anastasio, 2017). Specifically, the MS-model takes the form of a recurrent network in order to use a more efficient learning procedure to train its more extensive representations of neurobiological interactions, and to conform to a larger set of experimental observations.
Our approach is to model the interactions within and between the monoaminergic neurotransmitter systems and the stress hormone system as a network of nonlinear elements, and to train the model using machine learning. Both the structure of the model and the data on which it is trained are based on experimental observations as described in the literature. The model is of necessity abstract, but the form and behavior of the fully trained model are consistent with a broad range of findings on the monoaminergic-transmitter and stress-hormone systems. The fully trained model can be used to predict the acute effects of drug or drug-hormone combinations. Through systematic adjustments of model parameters that represent known, adaptable elements of the monoaminergic-transmitter and stress-hormone systems (referred to as transmitter system components, TSCs), the model can also be used to predict the range of responses of a population of patients to chronic administration of drug or drughormone combinations.
Acute administration of substances (such as drugs or hormones) has been observed to alter neuronal activity levels, and chronic (days to weeks) substance exposure can lead to adaptive changes in neurons that move their activity levels back toward their original levels (Blier and De Montigny, 1987;Turrigiano, 1999;2008). We simulated neuroadaptive changes by allowing a subset of TSCs (mainly proteins such as neurotransmitter or neurohormone receptors or transporters) to adjust their strengths (corresponding to factors such as expression levels, sensitivities, and synaptic locations) incrementally up or down. TSC-strength configurations that restored the activities of DR, LC, VTA, and PVN back toward normative baselines were referred to as "adapted. " We found that many different TSC-strength configurations could produce adaptation to chronic administration of any drug or drug combination, and each adapted configuration had its own, unique pattern of neurotransmitter and hormone levels. This computational finding provides a possible explanation for the clinical finding that no single antidepressant drug or combination alleviates depressive symptoms in all patients (Barowsky and Schwartz, 2006;Turner et al., 2008;Zhou et al., 2015;Cipriani et al., 2016;Friedrich, 2017). By revealing that many different TSC-strength configurations, each associated with its own monoamine levels, could be equally well adapted to any chronic antidepressant treatment, the MS-model identifies a key challenge to effective antidepressant drug design. By enumerating a very large set of possible neuroadaptive configurations to chronic drug or hormone administration, and by reporting the associated monoaminergic neurotransmitter and cortisol levels, our computational model can be used to address this challenge and to identify drug and/or hormone combinations that potentially could be therapeutic for a higher proportion of patients, or patients with specific subtypes of depression requiring elevations in specific monoamines.

Model Formalism Overview
The MS-model takes the form of a recurrent network of units that all have the same sigmoidal (S-shaped, nonlinear but differentiable) activation function. This powerful computational formalism is a Turing-equivalent universal approximator that can be efficiently trained, via machine learning, to perform a broad range of desired, dynamic input-output transformations (Siegelmann and Sontag, 1991). When used to model neural systems, the units in a recurrent network usually represent single neurons, or the average activity of the neurons in the same brain region (Anastasio, 2010). In modeling and analysis of the transmitter and hormonal systems that mediate mood, however, multiple levels of organization must be taken into account that include not only brain regions and neurons but also proteins and small molecules.
In the MS-model, individual units represent the brain regions, transmitters, receptors, transporters, enzymes, precursors, metabolites, and hormones that are involved in the pathologies of anxiety and depression. Some individual units represented whole brain regions that are central to these pathologies. The monoaminergic neurotransmitter-producing regions (DR, LC, and VTA) were represented as single units because the majority of antidepressant drugs target these regions directly (Koenig and Thase, 2009). The HPA axis regions (PVN, pituitary gland, and adrenal gland) were represented as single units in order to incorporate the effects of the stress response in our analysis. The amygdala, prefrontal cortex (PFC), and hippocampus were represented because these regions are implicated in regulation of the HPA axis, and because they are also implicated in the antidepressant response itself through their involvement in cognitive control (Albert et al.,2014;Dinan, 1996;Malagie et al., 1996;Godlewska et al., 2012).
Single units also represented key neurotransmitters (e.g. 5HT, NE, and DA) and hormones (e.g. cortisol (CORT) and oxytocin (Oxt)), and many of the key transmitter receptors, transporters, and enzymes that are transmitter-system components (TSCs). See Supplemental Figure 1 for a diagram of the full MS-model. A comprehensive literature search guided model-structure design in order to ensure that model architecture conformed to known neurobiological interactions.
Representing neurobiological entities (brain regions, transmitters, receptors, etc.) as single units enabled the model to represent the level (of activity, expression, concentration, etc.) of that entity in the whole brain, or in specific brain regions as appropriate. The 5HT unit, for example, represented the brain serotonin level, which is the major neurobiological endpoint for antidepressant action. The weights of the connections between the units also had specific identities in the model. For example, the effectiveness of the 5HT autoreceptor in inhibiting DR neurons was represented in the model by the absolute value of the inhibitory weight of the 5HT1A receptor (5HT1AR) unit onto the DR unit.
Relative to the overall number of network weights, a very small number of weights were of great significance because they represented TSCs whose efficacies, or strengths (expression levels, sensitivities, concentrations, synaptic locations, etc.), are known to adapt under chronic stress or chronic antidepressant administration. All of the weights in a trained network represent the normative strengths of influences of specific neurobiological entities on each other. The weights of the adaptable TSCs specifically were further adjusted to analyze the possible modes of adaptation to chronic stress, drugs, or hormones in the model.

Model Structure and Function
In describing neural networks it is necessary to distinguish inputs/outputs to/from the network overall, and inputs/outputs to/from individual units in the network. Units are categorized as input, output, or "hidden units. " There were 40 input units that provided input to the overall network. They take on assigned values and project to other units but do not receive connections from other units. Hidden and output units receive connections from input units and from each other. There were 23 output units that provided the output of the overall network. Output units are distinguished from hidden units in having targets (desired outputs). Hidden units have neither assigned nor desired values. The network consisted of 102 total units (input, hidden, and output).
Individually, the activity of each unit is a function of its inputs from the other units (this is not true for input units that do not receive inputs from other units). Each sending unit provides inputs to receiving units that are equal to the product of the sending unit's activity level and the value of the weight (positive or negative) of the connection from the sending unit to the receiving unit. The net input to a receiving unit is the sum of the weighted inputs from all its sending units. The activity level of a unit is the value of its net input (i.e. weighted input sum) after it has passed through the sigmoidal activation function, which "squashes" the net input sigmoidally in the range between 0 and 1. The value of the sigmoidal squashing function for net input 0 is 0.50 (the midpoint of the 0-1 range).
The connection weights from the input, hidden, and output units to the hidden and output units are organized into the network weight matrix. The weight matrix consisted of 61 rows and 102 columns for a total of 6222 weights. Any input, hidden, or output unit could project to any hidden or output unit, but various distinctions between specific classes of connection weights were made in order to satisfy modeling goals.
The subject of the MS-model is the behavior of the monoaminergic transmitter systems and the stress hormone system (i.e. the monoamine-stress system) as they relate to the pathophysiology of depression. This behavior involves heavy interaction between the three monoaminergic transmitter systems and the stress hormone system, but also involves interactions between those systems and many other systems throughout the nervous system. The MS-model is therefore focused on the interactions between the monoaminergic transmitter and stress hormone systems, but also represents many other interactions that are included to ensure that model behavior is in line with a large set of empirical observations on the monoamine-stress system, as described in the literature.
In order to foreground monoamine-stress system interactions in the model, three classes of connection weights are distinguished: canonical, structure, and non-structure. Canonical weights are the weights of the known connections between units representing key components of the monoamine-stress system, such as the weight from the DR to 5HT, representing the effectiveness of the DR in producing 5HT. There were 23 canonical weights (see Supplemental Table 2). Structure weights are the weights of the connections between units representing neurobiological entities that are also known to interact empirically but are not canonical weights, such as the weight from DR to galanin, representing the effectiveness of the DR in co-releasing galanin. There were 305 total structure weights. Non-structure weights denote all other connection weights; they may or may not represent as yet unidentified interactions that actually do occur in neurobiology. There were 5,894 non-structure weights in total. The three classes of weights are treated differently during model training (see next subsection).
Two data structures were needed in order to construct and train the MS-model: a structure matrix and a truth table. Both the structure matrix and the truth table were compiled via a comprehensive literature search. To include as much of the available experimental data as possible, the search compiled findings from multiple groups using a broad range of experimental methods in humans and other animals. Due to the heterogeneity in methods and test subjects, values in the structure matrix and truth table do not reflect precise measurements (e.g. mg/kg of a compound in a rat versus a human) but were quantized into a few discrete values. This course-graining allowed the model to represent much of the relevant phenomenology but in an approximate way.
The structure matrix is a two-dimensional matrix, coextensive with the network weight matrix, which specifies which neurobiological entities in the model are known to interact with which others, and their valance (positive or negative) if also known (see Supplementary Material S1: Details on Structure Connections). The structure matrix designates all structure and canonical weights in the model using non-zero integers (+/− 1 for structure, +/− 2 for canonical). Non-structure weights take value 0 in the structure matrix. Figure 1 shows a highly simplified version of the structure of the model using only one of the three monoaminergic neurotransmitter systems in its representation. The full model diagram, which incorporates the connections within and between all three of the monoaminergic neurotransmitter systems and many other neurobiological entities, can be viewed in Supplemental Figure 1: Complete Model Structure Diagram.
The truth table is an array of input/desired-output training patterns that specifies how specific experimental manipulations, which are represented as patterns of network inputs, are known to affect specific neurobiological endpoints, which are represented as patterns of desired network outputs. Each row of the truth table represents the statistically significant results of one or more actual experiments. When more than one finding was available on a particular input-output relationship, a consensus was reached using the available data (see Supplemental Material S2: Truth-table Justification for a summary of the experiments included in the truth table and corresponding references).
Inputs in the truth table are either present or absent (1 or 0) and outputs are assigned discrete levels between 0.30 (maximal decrease) and 0.70 (maximal increase) with reference to a baseline of 0.50. An output could either decrease maximally (0.30), decrease moderately (0.40), exhibit no change from baseline (0.50), increase moderately (0.60), or increase maximally (0.70) in response to an input or combination of inputs. Determinations of "moderate" versus "maximal" changes were made based on a comprehensive literature search of available qualitative and quantitative data. The discrete level assigned to a neurotransmitter as an output was taken as the whole-brain level of that neurotransmitter, if reported. If neurotransmitter level was reported in specific brain regions, then its level in the PFC was used preferentially over other brain regions. Otherwise, the discrete level assigned to a neurotransmitter was proportional to the average of its level in the regions in which it was measured.  Table).

Training the Model
The model was trained using an efficient, gradient-based machine learning method known generally as recurrent back-propagation. The specific algorithm we used is due to Piñeda and it trains dynamic, recurrent neural networks to produce desired steadystate output patterns given specific input patterns (Pineda, 1987) (see Supplemental Material S3: Details on Model Training). The Piñeda algorithm assumes that networks reach steady-states and tends to train networks to reach steady states.
Prior to training, the complete network weight matrix, including canonical, structure, and non-structure weights, is randomized. Training occurs over 1 × 10 6 training cycles October 2019 | Volume 10 | Article 1215 Frontiers in Pharmacology | www.frontiersin.org (presentation of input/desired-output patterns). Input/desiredoutput patterns are presented to the network in random order during training. The steady-state network response to each network input is found after 100 iterations of unit updating (each unit computes the value of the sigmoid for its net input on each iteration). The differences between the desired and actual outputs of the units in the model were used to compute an error measure. The training method then computes the change for each weight using the error measure, and updates network weights on each training cycle.
Weight changes were scaled by a learning rate term before being applied to a weight. The learning rates were set to 1 for the canonical and structure weights, and to 0.10 for the non-structure weights, in order to disadvantage the non-structure connections because they are not known for certain to be involved in the behavior being modeled. All weights had an upper-bound at absolute value 10. All weights had a lower bound of 0 except for the canonical weights, which had a lower bound of 1 to ensure they exerted an influence on overall network performance. For more details on model training, see Supplemental Material S3: Details on Model Training.
We developed a pruning procedure to prune networks in order to eliminate unneeded non-structure connections. This was intended both to minimize the number of non-structure connections and to improve generalizability of model behavior. Generalizability was assessed by training the model on all single-manipulation inputs (e.g. single drugs) and testing on all combination inputs (e.g. drug combinations). The truth table FIGURE 1 | Simplified schematic representation of the model. The Monoamine-Stress model takes the form of a recurrent neural network of nonlinear elements (units) that represent neurotransmitter-producing regions, enzymes, neurotransmitters, hormones, and receptors. Each unit type in the model is represented using a different shape in this highly simplified model diagram, in which only one or two of each unit type is shown. The full model diagram incorporates all three monoaminergic neurotransmitter-producing brain regions and the stress hormone system, and can be viewed in Supplemental Material. Neurotransmitter and hormone producing regions are represented as triangles, neurotransmitters and hormones are represented as circles, protein molecules are represented as rectangles, and inputs are represented as rounded rectangles. Connections between model units can be excitatory (arrow) or inhibitory (tee). DR, dorsal raphe; AG, adrenal gland; 5HT, serotonin; CORT, cortisol; 5HTT, serotonin transporter; 5HT1AR, 5HT1A receptor; and GCR, glucocorticoid receptor. included 41 single-input patterns and 25 combination-input patterns. Our pruning procedure eliminates 3981 non-structure connections (68% of all non-structure weights) on average and slightly improves generalization. The process of optimizing the pruning method is discussed in detail in Supplemental Material S4: Pruning Methods. All further networks were trained using the following procedure: Networks with the full weight matrix (i.e. all canonical, structure, and non-structure connection weights) were trained on the full truth table (all single and combination inputs). Nonstructure connections were pruned at the sensitivity cutoff optimized for generalizability, and the pruned networks were then retrained on the full truth table.

Adjusting Adaptable TSCs
Most antidepressants are administered chronically, so a model designed to represent the behavior of the monoaminestress system in the context of depression must represent the responses to chronic-drug as well as to acute-drug administration. The MS-model can reproduce the effects of acute antidepressant administration by virtue of its training on the truth table, which includes many input/desired-output patterns that describe the effects of acute drugs. Because the monoamine-stress system is known to adapt under conditions of chronic drug (or chronic stress), the model must also account for neuroadaptation. Neuroadaptation is the process by which the overall activity levels (specifically the resting or spontaneous activity levels) of neurons in key brain structures are brought back to their normative baselines under conditions of chronic perturbation, such as chronic administration of a drug (see Introduction).
The four key brain regions in the MS-model are the four canonical neural structures: DR, LC, VTA, and PVN. These canonical brain structures are represented as individual units in the model. Neuroadaptation is known to occur through changes (i.e. adaptations) in the expression, sensitivity, efficacy, cellular or synaptic location, or other aspects associated with proteins (receptors, transporters, ion channels, etc.) that influence the activity levels of neurons in key brain regions (see Introduction). The 10 adjustable TSCs in the MS-model correspond to 10 actual TSC proteins that are known empirically to undergo adaptive changes under conditions of chronic manipulations within the purview of the model truth-table (see Supplemental Table 3: Adjustable TSCs). These 10 real TSCs are also known to play critical roles in the interactions between the three monoaminergic-transmitter systems and the stress-hormone system.
The strengths of cell-type specific TSCs are represented as individual weights in the MS-model. For example, the best known adjustable TSC in the context of antidepressant neurobiology is the 5HT1AR on DR neurons, which has been observed to "desensitize" (i.e. to become less sensitive to 5HT and therefore less effective in inhibiting the DR) (Blier and De Montigny, 1987;El Mansari et al., 2005). The TSC corresponding to the 5HT1AR on DR neurons is represented in the model as the weight of the connection from the 5HT1AR unit to the DR unit.
TSCs are network connection weights, but changes in these weights in the context of neuroadaptation are fundamentally different from changes due to neural network training. In neural network training, all connection weights can be changed in order to bring network outputs closer to specific desired, target outputs. In neuroadaptation, only a small subset of weights may change (i.e. those that correspond to key TSC proteins that are known to adapt), and these changes are directed not toward achievement of specific target outputs but to produce a more general restoration of the overall activity levels of neurons in key brain structures. For this reason, neuroadaptation was produced in the model simply by incrementally changing (i.e. adjusting up or down) the strengths of the network weights corresponding to the 10 adjustable TSCs. Note that the weights corresponding to the strengths of the 10 adjustable TSCs constitute a small number of the network connection weights.
An adjusted network was one in which the ten weights corresponding to the adjustable TSCs are replaced with one or more adjusted values. Adjusted networks were considered adapted if their activation error was lower than their initial error, where initial error is the error due to acute administration of the drug or combination that is observed before any TSC adjustments have been made (see Results). The adjustment procedure was implemented "blindly" (i.e. in an unsupervised manner) so the adjustments could just as easily move the activities of the DR, LC, VTA, and PVN units away from normative activity as toward it, but adapted configurations were easily identified by the behavior of the adjusted networks.
We computed all possible configurations of the 10 adjustable TSC weights that were reachable by increasing or decreasing any single TSC weight by 0.50, within the predetermined TSC strength minimum of 0 and absolute maximum of 10, for a preset number of allowed adjustments that was the same for all weights. The "normative" strength of any network connection weight, including the weights corresponding to the 10 TSCs, is simply its value following neural network training. Due to the randomness inherent in neural network training (see previous subsection), equally well-trained networks can have very different network connection weights. This variability nicely corresponds to the natural variability in neurobiological properties that are known to occur between individuals (see Discussion).
To account for this inter-individual variability of neurobiological systems we studied adjustments of the normative weights in three representative networks, each trained from a different random initial weight matrix according to a different random order of input/desired-output presentations. We found all configurations of the 10 TSC weights reachable from each of the three representative, normative networks for all increments of 0.50 within the bounds of 0 and |10| up to a total of six adjustments, producing 382,747 total TSC-strength configurations over the three networks. Possible modes of adaptation of the MS-model were studied by analyzing this large set of adjusted configurations. Generation of the full set of TSC-weight configurations for a total of seven adjustments, which would have produced over 11 × 10 6 configurations, was not possible due to technical limitations (see subsection on Hardware Considerations).

Full-Range Individual-Weight Adjustment
We assume that all normal individuals will adapt to chronic antidepressant but we further assume that different individuals will adapt in different ways, and not all of them will achieve "therapeutic" levels of the monoamines, defined as the levels necessary to achieve remission of depressive symptoms. We set the therapeutic floor for 5HT to 0.70 because the baseline for 5HT was set to 0.50, and the desired 5HT output for acute SSRI was set to 0.60 (see Supplemental Material S2: Truth Table Justification). The therapeutic 5HT floor was set to 0.70 in order to double the acute increase of 0.10. This is consistent with findings that 5HT levels can increase to about 200% and 400% of baseline with acute or chronic SSRI administration, respectively (Ceglia et al., 2004). Elevations in NE and DA with chronic antidepressant administration have also been found to be associated with antidepressant response (Page et al., 2003;Yoshitake et al., 2004). For consistency, 0.70 was also used as the therapeutic floor for NE and DA because therapeutic levels for NE and DA with chronic antidepressant have not been identified conclusively. CORT levels have been found to decrease with chronic SSRI administration and antidepressant response, so the therapeutic ceiling for CORT was set to 0.70 to reflect a decrease from the acute-SSRI CORT level of 0.70 (Lenze et al., 2011;Ruhé et al., 2015).
It is possible that some adjustable TSCs contribute more to the attainment of therapeutic monoamine levels than others. We define "therapeuticity" as the ability of a TSC to contribute to the attainment of therapeutic monoamine levels through adjustments in its strength. We conducted full-range individualweight adjustment (FRIWA) analysis (a kind of sensitivity analysis) to gauge the therapeuticity of single TSCs. The starting configurations for FRIWA analysis were all of the configurations adapted to chronic SSRI that were also therapeutic, extracted from the exhaustive set of configurations of adjustments in all TSCs up to a total of six adjustments. The total number of adapted and therapeutic configurations over all three representative networks was 5598. Therapeutic configurations were those that increased 5HT above the 5HT therapeutic floor (> = 0.70) and decreased CORT below the therapeutic CORT ceiling (< = 0.70).
FRIWA analysis involved adjustment of the weight of a single adaptable TSC across its full range (0 to |10|; note that individual TSCs are either positive or negative), while the weights for the nine other adaptable TSCs remained frozen at their starting values, in each of the 5598 adapted and therapeutic starting configurations. FRIWA occurred in steps, where each individual weight adjustment (IWA) was an increment of |0.50|. Thus, FRIWA generated a set of 200 new configurations (20 IWA adjustments for each of 10 TSC weights) starting from each of the 5598 adapted and therapeutic configurations for a total of 980,193 new configurations. All configurations that were no longer adapted after a step of IWA were excluded from further analysis, leaving 493,564 adapted TSC-strength configurations. Of those, 285,635 configurations were designated "resistant, " because they remained therapeutic despite a step of IWA, while the remaining 207,929 configurations were designated "sensitive" because there were rendered non-therapeutic by a step of IWA.
The weights of each of the 10 adjustable TSCs in all of the post-IWA configurations were pooled over each representative network, and separated on the basis of resistance and sensitivity, excluding the weights that were manipulated by FRIWA. The mean weight of each adjustable TSC was computed for both the resistant and sensitive configurations of each network and compared (see Results). Further FRIWA analysis involved computing the pairwise correlations between all TSC weights over all the resistant, or over all the sensitive, configurations in each network. Only pairwise correlations that were statistically significant (P < 0.05) over all three representative networks would have been reported. Again for correlation analysis, the TSC weights that were manipulated by FRIWA were excluded.

Temporal-Logic Model-Checking
As described above, we studied the consequences of adaptation by generating a large set of adapted TSC-strength configurations (see subsection on Adjusting Adaptable TSCs). We also studied the process of adaptation by making allowed TSC weight adjustments (increments of 0.50 up or down, within bounds of 0 and |10|) in all possible sequences. Linear temporal logic (LTL) is a type of modal temporal logic that facilitates reasoning about sequences of discrete states evolving in time. LTL analysis enables the evaluation of temporally specified logical propositions such as whether a specific state is always maintained or eventually reached; whether a specific state pertains only until another state pertains; or whether a specific state always leads to another specific state.
Temporal-logic model-checking allowed us to determine such temporal relationships between TSC-strength configurations (i.e. states) of the MS model. We used LTL model checking to determine if specific degrees of neuroadaptation (e.g. a configuration in which a specific TSC has been adjusted down three times) always leads to an adapted and therapeutic configuration (for all possible sequences of adjustments proceeding from that configuration up to a total of six adjustments) (see Supplemental Material S5: Details on Temporal-logic Model-checking Procedure for details on model-checking statements).
In order to evaluate model-checks in the MS-model, it was necessary to define a set of logical predicates to be tested. The following predicates were used in the temporal-logic analysis: fht_high, 5HT is above the 5HT therapeutic floor (> = 0.70); cort_low, CORT is below the therapeutic CORT ceiling (< = 0.70); TSC_sens_gt_3, the adjustable TSC has sensitized by at least three steps; and TSC_desens_gt_3, the adjustable TSC has desensitized by at least three steps. Note that TSC in the last two predicates is a placeholder for any specific, adjustable TSC (e.g. 5HT1AR). Then we evaluated the following propositions for each of the 10 adjustable TSCs, where |-> and /\ are the LTL "LEADS TO" and "AND" operators, respectively: These are equivalent to the LTL propositions that if at any point in the trajectory of TSC strength adjustments, the TSC sensitizes or desensitizes by at least three increments, then 5HT is high and CORT is low at some subsequent point in time. We found that October 2019 | Volume 10 | Article 1215 Frontiers in Pharmacology | www.frontiersin.org both these propositions were false for all of the tested TSCs in all three representative networks (see Results).
Next, we evaluated the following more-easily-satisfiable propositions, where \/ and ~ are the LTL "OR" and "NOT" operators, respectively: These are equivalent to the LTL propositions that if at any point in the trajectory of TSC strength adjustments, the TSC sensitizes or desensitizes by at least three adjustments, then either 5HT is high and CORT is low at some subsequent point in time, or the TSC is no longer sensitized or desensitized. If one adjustable TSC was solely responsible for modulating the 5HT and CORT levels, then the model check for either of those two propositions for that specific, adjustable TSC should return True. Each LTL model check was carried out using each of the three representative networks for up to six adjustments, and only model-checking results that were consistent over all three networks are reported in Results.

Hardware Considerations
MATLAB was used for training, pruning, enumerating, and analyzing all adjustable-TSC configurations with up to six adjustments in TSC strength, generating all histograms and heatmaps, and FRIWA analysis. All MATLAB computational procedures were performed on an Intel Core 2 Duo CPU processor with 2, 2.33 GHz cores and 4 GB of RAM under the Windows 7 operating system, an Intel-inside CORE i7 processor with 2, 2.69 GHz cores and 8 GB of RAM under the Windows 8 operating system, and an Intel Core i7 processor with 4, 4.00 GHz cores and 32 GB of RAM under the Windows 10 operating system. Training a network (100 time steps, 1 × 10 6 iterations, 66 training patterns) took between 20 and 25 min on these machines. Exhaustively computing full sets of adjustable TSC configurations to degree 6 took about 1 minute. Computational overhead (memory limitations) prevented computing the full set to degree 7. Checking the set of TSC-strength configurations for neuroadaptation to chronic drug or hormone combinations in each of the three representative networks took about 8 min per network using MATLAB.
Python was used for enumeration of TSC-strength configurations and LTL analysis. Python search and LTL analysis of TSC-strength configurations having 127,582 states took 34 seconds on one quad-core Intel Core i5 processor. For subsequent model checks where steady-state values had been cached, model checks took an average of 3.10 seconds. Python enumeration of TSCstrength configurations was limited to degree 6 for consistency with the MATLAB analysis and because the results were unlikely to be different for degree 7 than for degree 6 (see Results).

Agreement Between Actual and Desired Outputs
Networks were trained, pruned, and then retrained as described in Methods. The results of training can be viewed in Figure 2.
Each plot in Figure 2 shows all of the desired and actual outputs for one brain-region, transmitter, or hormone output unit (DR, LC, VTA, PVN, 5HT, NE, DA, or CORT). Each actual output response and desired output response is plotted as a solid line and a dashed line, respectively. Correspondence between steady-state actual responses and desired responses is nearly exact. The results of training are shown in Figure 2 for only eight model units, however, the degree of correspondence between actual and desired output responses shown in this figure is representative of the remaining trained output responses. The average error of these units over all of the trained-prunedretrained networks was very low (6.39 × 10 −5 ), indicating close correspondence between model output responses to acute inputs described in the neurobiological literature. The heatmaps in Supplementary Figure 2 illustrate the extent of model agreement for all 23 output units over the whole training set both before ( Supplementary Figure 2A-C) and after pruning (Supplementary Figure 2D-F).

Enumeration of Adjustable TSC-Strength Configurations
The baseline activity levels of key model units are represented as dashed blue lines in Figure 3. The baseline activity is simply the activity of units in a trained network when the inputs to the network are all zero. All units have positive (nonzero) baseline activities because the squashing function, which determines each unit's output as a function of its net weighted input, produces a "spontaneous" output activity of 0.50 for a net input of 0. Note that the squashing function bounds unit activation between 0 and 1 (see Methods).
The units in the network influence each other's activity through their weighted interconnections. Because the unpruned weights of the connections between the units are nonzero (positive or negative), the baseline activities of the units can depart substantially from 0.50. The baseline activity of the units is basically the "response" of the network to 0 input. Following an initial transient, all unit "responses" to 0 input (baseline responses) settle into a stable activity pattern within 25 time steps that is maintained for the duration of the response. The baseline activity levels of the units in a trained network are their "normative" activity levels.
Simulated acute administration of an SSRI (an SSRI input of level 1) alters the responses of the canonical units in the model, as shown by the orange lines in Figure 3. Note that the responses of some of the units to acute SSRI were trained directly because these responses had been reported in the literature and were included in the truth table. The acute SSRI responses of the other canonical units are essentially estimated by the model, on the basis of its structure and on its training on the truth table overall. In general, unit activities can deviate substantially from their normative baselines due to acute SSRI. From the neuroadaptive standpoint, any deviation from normative baseline is considered an error that should be corrected by the adaptive process. We define "adaptation error" as the sum of the absolute differences from baseline of the activities of the DR, LC, VTA, and PVN units, because these are the canonical units in the MS-model. October 2019 | Volume 10 | Article 1215 Frontiers in Pharmacology | www.frontiersin.org FIGURE 2 | Agreement between desired (i.e., target) and actual outputs after pruning and re-training. All of the desired and actual outputs, collected over all input/ desired-output patterns, are represented in a single plot for each of the brain region and transmitter or hormone output units (DR, LC, VTA, PVN, 5HT, NE, DA, and CORT). Each output response is plotted as a solid line and each target (i.e. desired) output is plotted as a dashed line. Each of the outputs reach a steady-state value within 25 time steps. The RMS error of this network is 5.10 × 10 −5 . Note that the solid line (actual output) is superimposed on the dashed line (desired or target output), illustrating the accuracy of the training method.
FIGURE 3 | Model element activities in the baseline (no-drug) condition, acute (no-adaptation) SSRI condition, and chronic (adaptation) SSRI condition. Each plot corresponds to a different model unit as labeled. The blue dotted line in each plot shows the baseline activity level of a unit in the normal (no-drug) baseline condition. The red line in each plot shows the unit activity in the acute (no-adaptation) SSRI condition. Note that acute SSRI changes the activity levels of all of the units. The yellow line in each plot shows the adapted activity of a unit in an example, adapted configuration with chronic SSRI. Note that the adapted DR and PVN unit responses return closer to baseline, and the adapted 5HT and CORT responses increase and decrease, respectively. DR, dorsal raphe; LC, locus coeruleus; VTA, ventral tegmental area; PVN, paraventricular nucleus of the hypothalamus; 5HT, serotonin; NE, norepinephrine; DA, dopamine; and CORT, cortisol.
Chronic (as opposed to acute) administration of a drug (or combination) was simulated simply by keeping the input unit corresponding to that drug (or drugs, for a combination) at 1. We further defined "initial error" as the adaptation error of a trained network subjected to chronic administration of a drug (or combination) in the absence of any adjustments of the strengths of the weights representing the adjustable TSCs. An "adapted network" was any network that, due to one or more adjustments in TSC weights, had adaptation error lower than initial error. Model responses to acute SSRI, and adaptation to chronic SSRI, were of particular interest because SSRI is the drug class most often prescribed for the treatment of anxiety and depression. The responses of the canonical units in an example network adapted to chronic SSRI are represented in Figure 3 as yellow lines. This figure shows that adaptation to an SSRI can return the responses of DR and PVN back toward normal while also increasing 5HT levels and decreasing CORT levels, which is in agreement with experimentally and clinically observed chronic SSRI effects (Rutter et al., 1994;Ceglia et al., 2004;Nikisch et al., 2005;Rush et al., 2006).
Many other adapted networks, however, did not also have this pattern of adapted behavior. Transmitter and hormone responses to chronic drug administration, therefore, must be evaluated in many different TSC-strength configurations, derived from adjustments in the TSC weights of several different initial networks. As described above, our results are based on three representative networks (see also Methods).
We analyzed all of the configurations, starting from each of the three different, representative networks that could be generated by up to six adjustments at an increment of 0.50 (the number of adjustments was limited to six due to computational overhead, see Methods). The distributions of the three monoaminergic transmitter and CORT levels corresponding to all the adapted configurations of the three networks to chronic SSRI are shown as histograms in Figure 4. Each histogram in Figure 4 shows the numbers of configurations adapted to chronic SSRI that had levels of 5HT, NE, DA, and CORT falling within various bins as indicated. The first three columns show the histograms for the three representative networks separately, while the fourth column shows the histograms for the three representative networks combined. In this histogram figure and in all subsequent histogram figures the baseline, average, and therapeutic levels of each transmitter or hormone will be represented with a blue, green, or magenta line, respectively.
A central goal of this study was to use the trained network models to evaluate the possibility that combinations of antidepressant drugs, or combinations of drugs and hormones, could be more therapeutically effective that single SSRIs. Configurations adapted to chronic administration of select drugs and hormones in combination with SSRI were analyzed. The results of some specific drug and hormone combinations are shown in Figures 5 and 6. The rows of these two sets of histograms (Figures 5 and 6) show the monoamine and CORT distributions of all configurations starting from all three representative networks that were adapted to SSRI alone, to SSRI paired with another drug or hormone, and SSRI combined with the two other drugs or a drug and a hormone. Each column of these two FIGURE 4 | Histograms showing numbers of configurations adapted to SSRI expressing different levels of monoamines or CORT. Results are shown for three networks individually or pooled. Networks were adapted to chronic SSRI. The bin width of these and all other histograms was set to 0.03. The blue, green, or magenta vertical line in each plot is located at the baseline level (0.50), the average level, or the therapeutic cutoff (0.70), respectively, for each neurotransmitter or hormone. This figure also demonstrates pooling of the adapted configurations from three networks. October 2019 | Volume 10 | Article 1215 Frontiers in Pharmacology | www.frontiersin.org sets of histograms (Figures 5 and 6) corresponds to a specific transmitter or hormone unit (5HT, NE, DA or CORT). Viewing the histograms down the columns shows how the distributions of monoamine and CORT levels are affected by the different drugs and combinations. These histograms illustrate the wide degree of variability in adapted neurotransmitter and hormone levels that may be present in real, neurobiological systems with chronic drug or hormone administration. Figure 5A shows the results of adaptation to SSRI alone, 5B shows SSRI paired with Asenapine (an antipsychotic drug), 5C shows SSRI paired with Oxytocin (a peptide neurohormone), and 5D shows SSRI combined with both Asenapine and Oxytocin. Figure 5B shows that combining SSRI with Asenapine not only increases the proportion of adapted states with therapeutically elevated (toward the right) 5HT, but also shifts the NE histogram to the right as well, increasing the proportion of adapted states with elevated NE. Figure  5B also shows that the combination of SSRI and Asenapine can decrease CORT levels in a higher proportion of adapted configurations below the therapeutic ceiling (toward the left). In Figure 5C, the combination of an SSRI and Oxytocin also shifts the 5HT distribution to the right and increases the proportion of adapted states with low CORT. The combination of all three factors (SSRI, Asenapine, and Oxytocin) in Figure  5D increases the proportion of adapted states with high monoamine levels for all three monoamines, and also increases the proportion of adapted states with low CORT beyond the levels observed with SSRI by itself. These findings suggest that augmentation of SSRI action with Asenapine, Oxytocin, or both can potentially be therapeutic in a larger proportion of depressed patients than an SSRI alone. Figure 6A shows the results of adaptation to SSRI alone, 6B shows SSRI paired with Bupropion (an atypical antidepressant), 6C shows SSRI paired with Olanzapine (an antipsychotic drug), and 6D shows SSRI combined with both Bupropion and Olanzapine. Figure 6B shows that adaptation to the combination of SSRI and Bupropion can increase the proportion of adapted states with high NE and DA levels over that observed with the SSRI by itself. It also shows that this combination decreases CORT levels in a higher proportion of adapted states than the SSRI by itself. Figure 6C shows that the combination of SSRI and Olanzapine can increase the proportion of adapted states with therapeutic 5HT over the SSRI by itself. The combination of all three (SSRI, Bupropion, and Olanzapine) in Figure 6D was found to shift all three of the monoamine (5HT, NE and DA) histograms to the right (toward high monoamine levels) and increase the proportion of adapted states that reduce CORT below its therapeutic ceiling.
In none of the histograms in Figures 5 and 6 are the monoamine or CORT levels in all of the adapted states in their therapeutic ranges. This is in agreement with the finding that no single antidepressant drug, hormone, or combination has been found to be therapeutic in all patients. Overall, these histograms illustrate how combinations of chronic drugs (or of drugs and hormones) can increase the proportion of adapted states with elevated monoamines and reduced CORT levels (see (D) shows that combining an SSRI with both Oxt and Asn further increases the proportion of high monoamine and low CORT states. These histograms suggest that combining an SSRI with either Oxt, Asn, or both may be therapeutic for a greater proportion of patients than an SSRI administered alone. October 2019 | Volume 10 | Article 1215 Frontiers in Pharmacology | www.frontiersin.org Discussion). The shifts in neuroadapted monoamine levels with chronic administration of the drug and hormone combinations in Figures 5 and 6 correspond with experimental findings that these combinations may result in elevations in brain monoamines in rodent models (Shrier et al., 2000;Björkholm et al., 2014;Amini-Khoei et al., 2017;Thomas et al., 2017). These findings represent the model prediction that the combinations of drugs and hormones in Figures 5 and 6 may be therapeutic in a larger proportion of depressed patients than SSRIs alone, and are in agreement with clinical observations (Leuchter et al., 2008;Zhou et al., 2015;Pilkinton et al., 2016).

Predicting the Effects of Drug and Hormone Combinations on Monoamine Levels
The model was then used to evaluate the configurations adapted to chronic administration of 60 drug or drug/hormone pairs and triples. The 27 pairs consisted of SSRI paired with each of the 27 drugs and hormones (excluding SSRI itself) that were used as individual (single) inputs in model training. The 33 triples consisted of the subset of those pairs that have been used to augment SSRI action clinically or (in one case) experimentally, combined with a third substance that was either Oxytocin (a hormone), Antalarmin (a CRF1 receptor antagonist), or Olanzapine (an antipsychotic). These three substances were selected due to recent preliminary evidence suggesting that Oxytocin, Antalarmin, and Olanzapine could potentially have antidepressant effects either by themselves or in combination with another antidepressant (Björkholm et al., 2015;Scantamburlo et al., 2015;Amini-Khoei et al., 2017;Thomas et al., 2017).
The monoamine levels for the three networks adapted to each of the drug or hormone pairs and triples were compiled and the average monoamine levels were computed. The average levels of each of the three monoamines for each chronic drug or hormone combination are shown as rows in the heatmap in The excess monoamine reference vector represents a tripling of the acute level and is included in order to identify drug pairs and triples that may elevate the monoamines high enough to produce unwanted side effects (Shrier et al., 2000;Boyer and Shannon, 2005). All adapted monoamine vectors [5HT NE DA] that had any monoamine elevated above 0.80 were ordered by their vector distance from the excess monoamine vector; all remaining adapted monoamine vectors were ordered by their vector distance from the therapeutic reference vector.
The majority of the drug/hormone combinations tested in the model have been tested neither in animals nor in humans, but  D) shows that combining an SSRI with both Bup and Olan further increases the proportion of high monoamine and low CORT states. The histograms in (D) illustrate that the combination of SSRI+Bup+Olan shifts the monoamine and CORT histograms, especially those of NE and DA, toward more therapeutic states, suggesting that this combination can be therapeutic for a greater proportion of patients than an SSRI administered alone. a few have been examined in clinical studies and the results are in general agreement with the model. SSRI alone was found to modestly elevate monoamine levels chronically. The combination of an SSRI and Reserpine was found to be at the bottom (low adapted monoamine vector) of this figure, as expected due to the monoamine-depleting effect observed with Reserpine (Cooper et al., 1994). The combinations that include both an SSRI and Olanzapine are at the upper (more therapeutic) range of Figure 7, which is in line with the clinical finding that Olanzapine can be used to effectively augment SSRI action (Björkholm et al., 2015;Zhou et al., 2015;Thomas et al., 2017). The combinations of SSRI/Bupropion and SSRI/Quetiapine, have both shown to be effective in depressed patients and are also at the upper range of Figure 7 (Adson et al., 2004;Baune et al., 2007;Papakostas et al.,  2007; Leuchter et al., 2008;Cutler et al., 2009). The heatmap in Figure 7 illustrates how the MS-model could be used to screen chronic drug/hormone combinations for potential efficacy as monoamine-elevating treatments.

Full-Range Individual-Weight Adjustment (FRIWA) to Evaluate the Therapeuticity of Each Adjustable Weight
Though we identified and examined adjustments in 10 adjustable TSCs, it is possible that single TSCs can alone determine the therapeutic state of the MS-model, regardless of the values of the weights representing the other adaptable TSCs. We define "therapeuticity" as the ability of a biological factor to alter the properties of a biological system in a therapeutic direction. Specifically here, therapeuticity is the ability of a TSC to alter three of the properties of interest in the MS-model: the level of adaptation, the level of 5HT, and the level of CORT. Identification of the TSCs that could alone determine therapeutic state could enhance antidepressant drug design by identifying specific receptors or transporters that could be targeted in single-drug therapies.
In order to determine if therapeutic effects were specifically dependent on single adjustable TSCs, an analysis was conducted that evaluated the contribution of each individual, adjustable TSC to therapeutic adaptation with chronic SSRI administration. Therapeutic states were defined as adapted states that increased 5HT up to or above the therapeutic 5HT floor (> = 0.70) and decreased CORT down to or below the therapeutic CORT ceiling (< = 0.70). All configurations, generated from adjustments of all adjustable TSC weights from one to six adjustments, which were adapted to chronic SSRI and also therapeutic, were pooled for each of the three representative networks for full-range individual-weight adjustment (FRIWA) analysis.
In FRIWA analysis, for each adapted and therapeutic configuration, one adjustable TSC weight was adjusted across its full range (0 to absolute value 10), while the nine other TSC weights remained frozen. FRIWA occurred in 20 IWA steps of 0.50 each (see Methods). All states that were no longer adapted after an IWA step were excluded from the analysis. The stilladapted configurations that were resistant (remained therapeutic after an IWA step) or sensitive (became non-therapeutic after an IWA step) were then identified. The weights for each TSC were compiled for all resistant or all sensitive post-FRIWA configurations for each representative network. The average values of the weights for each TSC over either the resistant or the sensitive configurations were then computed for each network, but the TSC weights that were adjusted using FRIWA were removed because their values were manipulated. Figure 8 shows the average strengths of each adjustable TSC weight over all resistant or all sensitive configurations as an asterisk or diamond, respectively. Each adjustable TSC weight in Figure 8 has three asterisk-diamond pairs in three different colors (red, blue, and green), illustrating the results from each of the three representative networks on a single plot. This figure   FIGURE 8 | Comparison of average adjustable TSC strengths between resistant and sensitive configurations in all representative networks. Every configuration adapted to chronic SSRI that was also therapeutic (high 5HT and low CORT) at degree 6 was assessed for resistance to adjustments of each of the 10 adjustable TSCs. All TSC strength adjustments that resulted in configurations that were no longer adapted were excluded. Configurations that remained adapted and therapeutic following weight adjustments were determined to be "resistant" and adapted configurations that were no longer therapeutic following TSC adjustments were determined to be "sensitive." The average strength of each of the 10 adjustable TSCs in all of the configurations, excluding those in which that TSC itself was adjusted, was computed for both the resistant and sensitive configurations and plotted as asterisks or diamonds, respectively. The results for all three networks are represented on this single plot using three different colors (red, blue, or green) to distinguish between the mean adjustable TSC strengths of each network. Note that the average resistant and sensitive strengths for each adjustable TSC are very close in all three networks, illustrating that each individual TSC can provide a contribution to therapeutic resistance, but no single TSC by itself determines the therapeutic state. 5HT1AR, 5HT1A receptor; AR2, α-2 adrenergic receptor; GCR, glucocorticoid receptor; DAT, DA transporter; NET, NE transporter; 5HTT, 5HT transporter; CRF1R, CRF1 receptor; DR, dorsal raphe; LC, locus coeruleus; VTA, ventral tegmental area; PVN, paraventricular nucleus of the hypothalamus; Pit, pituitary gland; Adr, adrenal gland.
illustrates closeness in the averages of all the TSC strengths between the resistant and sensitive configurations in all three representative networks, demonstrating that no single TSC mediates the therapeutic state by itself. The results of this analysis suggest that all adjustable TSCs contribute to the attainment of the therapeutic state and that no single TSC is determinative of it. This model result suggests that there may not exist a single, individual TSC that can be targeted to alleviate depressive symptoms in all patients (see Discussion).

Pairwise Correlations Between Adaptable TSC Strengths
The FRIWA analysis considered one adjustable TSC at a time. It does not rule out the possibility that therapeuticity may be a property of correlations between pairs of TSCs rather than of single TSCs by themselves. To evaluate this possibility, pairwise correlations between all TSC weights over either all resistant or all sensitive configurations were computed for each representative network, but again, the TSC weights that were adjusted using FRIWA were removed because their values were manipulated. Even at the permissive significance level of P = 0.05, the analysis found no significant pairwise correlations between adjustable TSCs that were consistent over all three representative networks in either the resistant or the sensitive configurations. This result suggests that there are no two TSCs that can be targeted together to ensure therapeutic neuroadaptation in all patients. It exemplifies the challenge faced by investigators designing chronic drug and hormone combinations for relief of depressive symptomology.

Temporal-Logic Model-Checking
As described in previous subsections, the FRIWA and correlation analyses examined a large space of configurations adapted to an SSRI and determined that no single TSC, nor pairs of TSCs, mediates therapeutic resistance by themselves. Neuroadaptation to chronic antidepressant, where each TSC change is incremental, can also be examined as a process to determine if specific states or degrees of neuroadaptation must be reached prior to arriving at therapeutic configurations. Linear temporal-logic (LTL) analysis can be used to elaborate all pathways of sequential TSC-strength adjustments in order to evaluate possible temporal relationships in neuroadaptation (see Methods). This mode of analysis can be used to determine whether certain numbers of TSC-strength adjustments, once attained, will always lead to an adapted and therapeutic state.
The antecedent of the LTL propositions we analyzed specified that a specific TSC had been adjusted three times, either up or down (see Methods for details). The antecedent of three adjustments was chosen because it was halfway between 0 and 6, the total number of adjustments to which we were limited for technical reasons. The consequent queried whether all subsequent sequences of up to three adjustments in any subset of the 10 TSCs would all lead to an adapted and therapeutic state (5HT > = 0.70 and CORT < = 0.70). Each LTL model-check was carried out using each of the three networks for six total adjustments in all TSCs, and only those model-checking results that were consistent over all three networks are reported. All propositions returned false.
We next examined the more easily satisfiable propositions that once a specific TSC had been adjusted three times, either up or down, then all subsequent sequences of up to three adjustments in any subset of the 10 TSCs would all lead to an adapted and therapeutic state, or the TSC no longer maintains its degree of adjustment. This proposition allowed for the possibility that failure to maintain an adapted and therapeutic state occurred because the specific TSC had not maintained the specified level of adjustment. All of these more satisfiable propositions, however, again returned false, as shown in Table 2.
The LTL analysis shows that three adjustments of any single adjustable TSC weight up or down will not guarantee that a therapeutic state will be reached. All propositions returned false to degree 6, and it is unlikely that they would be true were the LTL analysis extended to greater degrees of adjustment because that would entail additional opportunities for TSC de-adjustment and overall de-adaptation. We did not attempt an exhaustive LTL analysis because the number of potentially relevant propositions that could be tested is simply too many. The result we generated, that three adjustments either up (sensitization) or down (desensitization) in no single TSC by itself can lead to states that are all therapeutic within three additional adjustments over the other nine TSCs or itself, supports the FRIWA analysis finding that adjustment of no single weight individually can determine the therapeutic state when the other TSCs can also adjust. These findings impose potential challenges to effective antidepressant drug design (see Discussion).

DISCUSSION
A consensus is forming around the general idea that multifactorial diseases should be treated with multidrug therapies (Perry et al., 2015;Xu et al., 2015;Anastasio, 2017). It has become common practice to treat psychiatric disorders with combinations of drugs, but the approach has been ad hoc and based mainly on clinical trial-and-error (Barowsky and Schwartz, 2006;Zhou et al., 2015). SSRI non-responders are often treated with combinations of two drugs , but systematic evaluation of the relative benefits of different combinations of two or more drugs in depression treatment has not occurred. The identification of novel, multidrug treatments for depression could proceed either through rational design or brute-force screening of drug combinations. The main challenge associated with rational design is the overwhelming complexity of the neurobiology of depression, which makes it extremely difficult to know a priori how any specific drug combination will work. The main challenge associated with brute-force screening is the sheer number of possible drug combinations, which grows geometrically with the number of drugs.
Our model addresses both challenges. It addresses rational design by computationally representing many aspects of the structure and function of the monoaminergic transmitter and stress hormone systems, whose interactions are central in depression neurobiology. It addresses combinatorial explosion by providing the means computationally to screen a large set of drug combinations, and to identify the most promising combinations for experimental evaluation. Addressing these dual challenges is crucial, because treatment with SSRI alone, which is currently the first-line treatment for depression, is less than 40% effective, and efficacy is not greatly increased when SSRI is combined with another drug (Turner et al., 2008;Zhou et al., 2015;Cipriani et al., 2016).
Our previous, initial model of the neurobiology of depression offered an explanation for the low SSRI efficacy rate. That model (the Monoamine or M-model) represented the monoaminergic transmitter systems (and some related transmitter systems) but did not represent the stress hormone system. It did, however, contain 11 adjustable TSCs. It also admitted of many different TSC-strength configurations, and many of those were adaptive in that adaptation error in the presence of simulated, chronic SSRI was lower than initial error. In the M-model, however, only 29% of the configurations adapted to chronic SSRI also elevated 5HT to therapeutic levels (Camacho and Anastasio, 2017). The conclusion was that the clinically observed efficacy of chronic SSRI is low (< 40%) because the real monoaminergic neurotransmitter system likewise has many ways to adapt to chronic SSRI, but only some of those ways also elevate 5HT to therapeutic levels.
Both the M-model and our current model (the MS-model) take the monoamine hypothesis as a starting point due to the finding that drugs that elevate monoamine levels are effective in treating depression in placebo-controlled trials (Turner et al., 2008). Monoamines may or may not be deficient in depressed patients (O'Reardon et al., 2004); however, depression relief is associated with drug-induced elevations in monoamine levels with chronic antidepressant administration (Haddjeri et al., 1998). Also, the antidepressant effects of rapid-acting treatment approaches (such as Ketamine and ECT) may involve the acute actions of these drugs on monoaminergic neurotransmission (Can et al., 2016;Pinna et al., 2018).
The leading alternative to the monoamine hypothesis posits that depression relief is associated with elevations in hippocampal brain-derived neurotrophic factor (BDNF) (Duman and Monteggia, 2006;Castrén and Rantamäki, 2010). Other hypotheses implicate neuroendocrine systems (Rubinow et al., 1998;Schmidt, 2005;Pariante and Lightman, 2008) or the relative activity levels of interacting neural systems (Mayberg et al., 2000;Albert et al., 2014). However, drugs that elevate monoamines also elevate hippocampal BDNF, influence neuropeptide transmitter systems, and alter the relative activity levels of interacting brain regions (Duman and Monteggia, 2006;Surget et al., 2008;Berger et al., 2010). We assume that therapeutically elevated monoamine levels may alleviate depressive symptomology through associations with previously described mechanisms of depression relief (Berton and Nestler, 2006;Stone et al., 2008;Andrews et al., 2015). The MS-model represents a significant advance over our previous, initial model (the M-model) in several respects. The MS-model represents the interactions of the three monoaminergic systems not only with each other but also with the stress hormone system. Model structure connections were selected on the basis of known interactions between model units, and the pruning procedure optimized generalizability of model predictions (see Supplemental Material S4: Pruning Methods). By adopting a recurrent nonlinear network formalism (see Supplemental Material S3: Details on Model Training), the MS-model represents many more relevant neurobiological details and conforms model behavior to a much larger set of empirical observations with improved agreement between observed and actual output values.
In the M-model, we analyzed the TSC-strength configurations of a single representative network, but since individuals can vary as well as their TSC-strength configurations, it is necessary to analyze more than one representative network. In the MS-model we therefore analyze the TSC-strength configurations of three representative networks. In the M-model, we analyzed all TSCstrength configurations reachable in three adjustment steps. In the MS-model we double the number of allowed adjustment steps to 6. This represents an increase in the number of configurations analyzed from 11,155 for the M-model (with 11 TSCs but only one representative instance) to 382,747 MS-model configurations (with 10 TSCs and three representative instances). In the current work we also extend our analytical repertoire to include graphical (distribution/histogram), sensitivity (FRIWA), correlation, and LTL analysis. The MS-model recapitulates and extends the key idea established in the M-model, that the brain has many ways to adapt to chronic drug administration but not all of those ways will achieve therapeutic goals. An innovation of the current work is the use of histograms to show how the monoamine and CORT levels are distributed over a large set of adapted TSC-strength configurations (Figures 4, 5, and 6). To the extent that the model reflects actual neurobiology, these histograms illustrate how we might expect a chronic drug or drug combination to affect a clinical population of people who differ in their pre-drug, starting configurations and also differ in their subsequent adaptation pathways. Specifically, they illustrate that we should not expect everyone in a depressed population being treated with chronic antidepressants to raise their monoamines, nor to lower their CORT, to therapeutic levels. Computing average monoamine levels over a whole distribution of computationally determined, adapted TSC-strength configurations (Figure 7), is a plausible way to predict the effects on monoamine levels of chronically administered antidepressant drugs or drug combinations.
The fact that only a subset of the many possible TSC-strength configurations is both adaptive and therapeutic raises the possibility that one, or only a few, of the adjustable TSCs are determinative of the therapeutic state. We used FRIWA analysis (our customized variant of sensitivity analysis) to explore this possibility for single TSCs in all three representative networks. We found that some configurations were "resistant" in that they remained adapted and therapeutic despite IWA applied to a TSC. Other configurations were "sensitive" in that they did not remain adapted and therapeutic despite IWA applied to a TSC. The average values of individual TSC strengths, however, were essentially the same in resistant versus sensitive configurations in all three networks (Figure 8), indicating that no TSC by itself is determinative of the therapeutic state.
We then used correlation analysis to explore the possibility that pairs of TSCs were determinative of the therapeutic state, and considered all pairs that were significantly correlated at the relatively permissive significance level of P = 0.05. We found that no pair of TSCs was significantly correlated in either the resistant or sensitive configurations over all three representative networks, showing that no pairs of TSCs are together determinative of the therapeutic state.
We then used LTL analysis to explore the possibility that a specific TSC, once it reached a certain level of sensitization or desensitization, could guarantee a therapeutic outcome. For technical reasons (see Methods), the number of steps of adjustment we allowed in the strength of each TSC was limited to 6, and our LTL analysis took this constraint into account. For each TSC, we tested the truth or falsehood of two LTL propositions essentially asserting that if a given TSC has reached a certain level of sensitization or desensitization, then the network is guaranteed to arrive at an adapted and therapeutic configuration no matter what the other TSCs do in the remaining adjustment steps. All of these propositions returned false in our LTL analysis ( Table 2). What all of these analyses show is that no single TSC, nor pair of TSCs, nor TSCs that have attained a certain degree of adjustment, can determine the therapeutic state. They demonstrate, more generally, that the properties of interest in the model (the level of adaptation, the levels of the three monoaminergic transmitters, and the level of CORT) are overdetermined by the 10 TSCs.
In its "overdeterminedness" (or degeneracy), the MS-model makes contact with other models depicting phenomena at both network and neuronal levels of neurobiological organization. In the vertebrate brain, sensory signals, motor commands, and information in general is represented not by single neurons but by networks of neurons. These representations are overdetermined because the pieces of information are few relative to the number of neurons in the network that represents them. In consequence, the same information can be distributed in many different ways over a network, and neurons in the same network can vary greatly in their response properties in what has been termed a non-uniform distributed representation (Anastasio, 1991;Anastasio, 2010).
Overdetermination of physiological properties by the values of multiple, relevant parameters has also been demonstrated for single neurons and computational models of thereof (Edelman and Gally, 2001;Golowasch et al., 2002;Bucher, 2005;reviewed in Marder and Taylor, 2011). The most well-known case in point was established by Eve Marder and colleagues in their model of the lateral pyloric neuron (LPN) in the lobster somatogastric ganglion (Taylor et al., 2009). In the invertebrate brain, information is sometimes represented not by neural networks but by single, very large (i.e. "giant") neurons that are identifiable from one animal to another. Experiments revealed that many different ion channels determine the electrophysiological properties of LPNs, and that the parameters that determine the function of ion channels of specific types can vary between LPNs from different decapods (Golowasch and Marder, 1992;Golowasch et al., 2002;MacLean et al., 2005;Schulz et al., 2006;. The Marder group developed a biophysically detailed model of the LPN and identified 17 parameters (i.e. the conductances and some related properties of 10 ion channel types) that plausibly could vary between individual LPNs. They generated ~6 x 10 5 different configurations of these 17 parameters and found, after evaluating them all computationally, that 1304 of those configurations endowed the LPN model with the same, realistic set of electrophysiological properties (number of spikes, spike frequency, spike duration, burst duration, etc.). Rather than grouping into specific sets of values, however, the Marder team found that LPN models having equally realistic electrophysiological properties could vary greatly in the values of their 17 defining, ion-channel parameters. By examining parameter value distributions and correlations, they found further that no single parameter nor pair of parameters was determinative of LPN model properties.
Analogously, we find that many different configurations of TSC-strengths can endow the MS-model with the property of adaptation to simulated, chronic SSRI, but we take this a step further and show that the adapted configurations can also differ in other properties, namely in their levels of the monoamines. Overall, these computational studies suggest that there is a large amount of degeneracy in neurobiological systems in general. The MS-model specifically implicates overdeterminedness as a consequence of the complexity inherent in the antidepressant response that poses a challenge to the identification of an antidepressant drug or hormone combination that is effective in all patients (Edelman and Gally, 2001). An approach that recognizes different subtypes of depression, and matches those with monoamine profiles predicted from computational models of the adapted antidepressant response, could provide an effective means of identifying promising drug combinations. The MS-model identifies drug and hormone combinations that could potentially be more therapeutic than single drugs for patients suffering from specific subtypes of depression. Because the therapeutic effects of chronic antidepressant use have been associated with elevations in monoamine levels and reductions in CORT levels, we were specifically interested in configurations that reproduced these experimentally and clinically observed changes (Nikisch et al., 2005;Lenze et al., 2011). Our model demonstrates that combining an SSRI with certain antipsychotics, atypical antidepressants, or hormones can increase the proportion of adapted states that are also associated with elevations in the monoamines and reductions in CORT.
The MS-model identifies Asenapine, Olanzapine, Bupropion, and Oxytocin as adjuncts to SSRI therapy that hold particular promise. Precedent for the use of these adjuncts has been established clinically. Physicians in practice frequently resort to the augmentation of SSRIs with antipsychotics such as Asenapine and Olanzapine in SSRI-monotherapy non-responders (Boulton et al., 2010;Han et al., 2013). The atypical antidepressant Bupropion has also been found to relieve depressive symptoms in SSRI non-responders and is widely used clinically, either by itself or in conjunction with SSRI Zisook et al., 2006;Leuchter et al., 2008). Intranasal administration of the hormone Oxytocin combined with chronic SSRI also can be more effective in treating depressed patients than SSRI by itself (Scantamburlo et al., 2015).
Results from the MS-model are consistent with these general findings in that pairing SSRI with either Asenapine, Olanzapine, Bupropion, or Oxytocin increases the number of adapted states with elevated monoamines and reduced CORT over that observed with SSRI alone (Figures 5 and 6). Furthermore, the triples composed of SSRI/Asenapine/Oxytocin or of SSRI/Bupropion/ Olanzapine further increase the proportion of adapted states associated with therapeutic changes in monoaminergic neurotransmitter and CORT levels in the MS-model ( Figures  5 and 6). More generally, the heatmap in Figure 7 shows that combinations that included an SSRI and an antipsychotic tended to be located at the upper range of monoamine levels (closer to the therapeutic reference vector [0.70 0.70 0.70]), and combinations that included Bupropion tended to elevate average NE and DA levels closer to the therapeutic reference vector. The average monoamine levels of combinations that included Oxytocin hormone were also located near the therapeutic reference vector.
The model shows a wide heterogeneity in effects on monoamine levels of different drug/hormone combinations, and this is of potential clinical significance as different subtypes of depression have previously been shown to respond best to elevations in the levels of the different monoamines (Parker, 2000;Malhi et al., 2005). Specifically, patients with non-melancholic, melancholic, or psychotic depression have been observed to respond best to antidepressants that elevate 5HT, NE, or DA, respectively (Parker et al., 1992;Guelfi et al., 1995, Malhi et al., 2002Malhi et al., 2005). Drug combinations that included Trazodone (an SNRI), Olanzapine (an antipsychotic), and GBR-12935 (a DAT antagonist) tended to raise 5HT, NE, and DA levels, respectively, in the model. The computational screen predicts potentially adverse as well as beneficial combinations. Notably, combinations that include both an SSRI and an MAOI were found to increase monoaminergic neurotransmitter levels well beyond that of other combinations. This MS-model result is in line with the clinical finding that combining these two drugs can result in toxic elevations of the monoamines (Remick and Froese, 1990;Shrier et al., 2000;Boyer and Shannon, 2005). Overall, the model identifies many possibly beneficial drug/hormone combinations that could be considered for clinical trial after extensive testing for safety and efficacy using animals.

CONCLUSION
Our overall MS-modeling strategy is both rational and brute force. It is rational in that the model incorporates, both in its structure and its function, many aspects of the known interactions between the monoamine and stress-hormone systems-two systems that are central to depression neurobiology. It is brute force in that it can be used to screen large numbers of drug combinations computationally, and to predict the distributions of monoamine levels that could be expected in a diverse patient population in which each individual can adapt to the chronic drug regimen in their own unique way. The MS-model could be used to identify drug/hormone combinations that, once verified in animal experiments, could be therapeutic for a higher proportion of patients than single-drugs by themselves.

AUTHOR CONTRIBUTIONS
MC performed all literature searches, wrote the MATLAB computer programs, performed computer simulations and analyses, collected data, and wrote the manuscript. WV developed LTL procedures, performed LTL analyses, and co-wrote the manuscript. TA designed the study, developed the formalism and training methods, wrote the initial computer programs, directed the research, and co-wrote the manuscript.

ACKNOWLEDGMENTS
We thank Ignacia Caviedes, Haven Comeaux, Kate Hamblen, Emily Hamm, Neena Joshi, Mahak Lalani, Diana Masolak, Cassandra Mora, and Katherine Zitello for their help in compiling the database of experimental facts on which the MS-model is based. We thank Allegra Domel for her contribution to the computational pruning analysis. We also thank Shivali Patel for her computational contribution to the construction of the full model-diagram.