Computational modeling of opioid-induced synaptic plasticity in hippocampus

According to a broad range of research, opioids consumption can lead to pathological memory formation. Experimental observations suggested that hippocampal glutamatergic synapses play an indispensable role in forming such a pathological memory. It has been suggested that memory formation at the synaptic level is developed through LTP induction. Here, we attempt to computationally indicate how morphine induces pathological LTP at hippocampal CA3-CA1 synapses. Then, based on simulations, we will suggest how one can prevent this type of pathological LTP. To this purpose, a detailed computational model is presented, which consists of one pyramidal neuron and one interneuron both from CA3, one CA1 pyramidal neuron, and one astrocyte. Based on experimental findings morphine affects the hippocampal neurons in three primary ways: 1) disinhibitory mechanism of interneurons in CA3, 2) enhancement of NMDARs current by μ Opioid Receptor (μOR) activation and 3) by attenuation of astrocytic glutamate reuptake ability. By utilizing these effects, simulations were implemented. Our results indicate that morphine can induce LTP by all aforementioned possible mechanisms. Based on our simulation results, attenuation of pathologic LTP achieved mainly by stimulation of astrocytic glutamate transporters, down-regulation of the astrocytic metabotropic glutamate receptors (mGlurs) or by applying NMDAR’s antagonist. Based on our observations, we suggest that astrocyte has a dominant role in forming addiction-related memories. This finding may help researchers in exploring drug actions for preventing relapse.


Introduction
Addiction is a chronic brain disease, in which addicted patients continually exhibit drug-seeking behavior despite its adverse consequences. In an addicted person, cellular and molecular adaptations in particular brain regions lead to behavioral abnormalities [1,2]. Drugs induce plasticity in the functions of neurons and synapses. These modifications can be seen as a form of "molecular or cellular memory" [3]. In the addiction process, abnormal memories are stored in specific regions of the brain, including the hippocampus, amygdala, prefrontal cortex, Ventral Tegmental Area, Nucleus Accumbens [2,3]. As a matter of fact, these abnormal memories a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 are responsible for addiction syndromes. Drug-induced changes in the functions of glutamatergic synapses have been considered to be a significant cause of drug-associated memories [2][3][4][5][6]. Contextual information of drug consumption is stored in the hippocampus and plays a vital role in drug-related reward property [6,7] and the relapse [8].
Opioids that are used for the treatment of pain are intensively addictive. Opioids similar to other drugs may have long-term effects on the affected neurons through Long-Term Potentiation (LTP) or Long-Term Depression (LTD) mechanisms which are thought to be the cellular basis of memory formation [1,2]. Many opioids, including morphine, affect the neurons through μ opioid receptors, which are G-protein coupled receptors (GPCRs). In the presence of morphine, μOR activation can lead to induction of LTP/LTD [9]. Therefore, assessing opioid-induced LTP/LTD could be the key to understanding the mechanism by which opioids induce addiction [2,10]. Some animal research indicates that morphine affects the hippocampal synapses by inducing LTP in the injected areas, and thereby may trigger the addiction process [11,12].
Opioid-induced induction of LTP in hippocampal neurons is via a disinhibitory mechanism which affects the hippocampal interneurons [13][14][15][16][17]. Opioids cause an inhibition of GABA release from hippocampal interneurons [16,17]. So, pyramidal neurons of CA3 receive less GABA, and inhibitory postsynaptic potentials (IPSPs) in CA3 neurons attenuate [18]. It seems that inhibition of interneurons by opioids is due to the enhancement of K + channel activity, or due to the reduction of Ca 2+ channel currents [19]. It has also been shown that inhibition of Ca 2+ channels by morphine requires G-proteins activation [20,21].
The goal of this paper is to present a biophysically-induced mathematical model at the synaptic level that can denote some features of the cellular addiction formation by opioids in hippocampus. Based on the experimental evidence, we analyze the effects of opioids on LTP at the synaptic level due to the activation of μORs using the proposed mathematical model.
In the first part of the paper, neurobiology of opioid signaling in synaptic plasticity and memory formation will be described. The second part introduces the modeling approach to describe the opioid signaling process at the synaptic level for the induction of pathological LTP. In the third part, simulation results are shown for normal and pathological conditions. Finally, discussion and conclusion of the results are presented.
Neurobiology of opioid signaling in synaptic plasticity and addiction memory formation channels of interneurons [19]. Inhibition of Ca 2+ channels in turn leads to less Ca 2+ influx into the interneurons. The consequence of less Ca 2+ entry is the suppression of GABA release from interneurons. Thus, GABAergic receptors of CA3 pyramidal neurons receive less GABA, which causes more neuronal excitation [18]. Further excitation of pyramidal neurons leads to more glutamate release, which means CA1 pyramidal neurons and astrocytes are going to be more excited than normal level. Furthermore, opioids decrease astrocytic glutamate trans-porter1 (GLT1) activation which is responsible for the glutamate uptake from the synaptic cleft [28,30]. As a result, reuptake of glutamate is reduced by astrocyte, and therefore more glutamate can be found in the synaptic cleft. Enhancement of glutamate in synaptic cleft conveys to more excitation of the postsynaptic CA1 neurons.
Opioids also affect the CA1 postsynaptic pyramidal neurons. In fact, it removes the blocking effect of Mg 2+ on N-methyl-D-aspartate (NMDA) receptors and thereby increases the fast glutamatergic transmission [14,31,35]. The outcome of this process is the augmentation of calcium entry into the postsynaptic neurons [36,37]. Enhancement of calcium concentration in the postsynaptic neurons results in phosphorylation of α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) receptors and Nitric Oxide (NO) synthesis by Ca 2+ Calmodulin-dependent protein kinase II (CAMKII) mechanisms [31,38]. Then, NO retrogradely diffuses to the presynaptic CA3 pyramidal neuron and increases the binding capability of presynaptic calcium sensors, which may cause more glutamate release and induction of LTP in the synapse [31,39,40].

Modeling approach
As mentioned earlier, the goal of this paper is to present a biophysically-driven mathematical model at the synaptic level in hippocampus in order to describe the cellular addiction formation by opioids. The main components of the proposed model are: (1) a bouton of a presynaptic CA3 pyramidal neuron, which can produce action potentials (APs) by receiving external stimulations and generates NO which can act as a retrograde signal that enhances the release of glutamate, (2) a GABAergic single compartment interneuron at the CA3 layer which inhibits the CA3 pyramidal neuron, (3) a spine of a postsynaptic CA1 pyramidal neuron which can produce Excitatory Postsynaptic Potentials (EPSPs) based on the amount of glutamate in the synaptic cleft; functions of NMDA and AMPA receptors of the postsynaptic neuron are also considered in this model, and (4) an astrocyte which monitors the synaptic transmission, regulates presynaptic glutamate release and is involved in the postsynaptic transmission by releasing gliotransmitters. The block diagram of different parts of the model is presented in Fig 2.

Presynaptic neurons
For modeling the presynaptic CA3 pyramidal cell and interneuron, we have utilized Hodgkin-Huxley (H-H) type formulation to generate APs in presynaptic neurons [41]: APs can be generated due to externally applied current (I app ) in the presence of sodium (I Na ), potassium (I K ) and leak (I L ) ionic currents (full description of f Pyr and f Int can be found in Appendix). Besides, we have adapted and added a GABAergic inhibitory current (I GABA A ) to the presynaptic pyramidal neuron similar to the formulation used by Zachariou et al. [42]: In Eq (3), g is the ratio of activated GABA A receptors and is described by [34,42]: Here, g 1 (V Int ) is assumed to be a sigmoid function which its half maximum voltage depends on activation of the interneuron calcium channels (CaCh). τ g = 1 ms represents time constant of GABA A receptor.
Hippocampal opioid receptors which are activated by morphine are mostly μORs. These receptors belong to GPCRs which are activated by morphine or its agonists. Because of opioids disinhibitory effect on the synaptic transmission in the hippocampus, it is supposed that presynaptic opioid receptors are on the hippocampal interneurons. Recently, Zachariou and colleagues proposed a model for hippocampal cannabinoid receptors (CBRs) [42]. This model is based on Bertram's minimal model of G-protein activation [34]. Since μORs and CBRs are both GPCRs, we have adapted this model to describe the activation of μORs in the hippocampus. The presence of morphine (Morph) results in activation of μORs (MOR) which modulates calcium channels activation (CaCh). Equations derived from [34,42] are modified as follows for opioid receptors: Maximum number of activated μORs is described by MOR 1 . This parameter depends on morphine concentration (Morph). k_ Shows reluctant to willing transition rate of calcium channel and depends on interneuron's voltage (V Int ), this parameter reflects the voltage-dependent dissociation of G βγ from the channel, k + is the transition rate of the willing to reluctant state for calcium channel and depends on morphine concentration, this parameter reflects the concentration of the activated G-proteins. It can be seen that morphine presence activates μORs which is described in Eqs (5) and (6). Consequently, calcium channels activation in interneuron changes in Eqs (7) and (8). Therefore, the conductance for inhibitory current varies in Eq (3) due to changes in Eq (4). Here, for the sake of simplicity, like Bertram et al., we assumed no dynamics for the neurotransmitter release from presynaptic interneuron to the pyramidal postsynaptic neuron [34]. Instead of that, we assumed that activation of VGCCs in interneuron affects directly the conductance of GABA A receptors in the presynaptic pyramidal neuron [34].
After APs are generated in presynaptic pyramidal neuron, glutamate will be released. Glutamate release depends on calcium concentration in the neuron. Calcium concentration in the presynaptic neuron can change through fast (c Fast ) and slow (c Slow ) dynamics which is described by Eq (9).
Fast calcium oscillations which occur due to Aps, are governed by: Here, I Ca represents calcium current through N-type channels, I PMCa denotes calcium current due to ATPase pump which pumps extra calcium out of the cell, and I PMleak is the leak current from extracellular space into Bouton (Details about f Fast can be found in the Appendix).
Slow calcium oscillations occur due to the activation of the extrasynaptic mGlurs. Since mGlurs activation has an indispensable role in synaptic plasticity governed by opioids, we have considered its role in the model. Activation of mGlurs is due to the binding of gliotransmitters. The presence of gliotransmitters lead to Inositol trisphosphate (IP3) production in an intracellular space which enhances calcium concentration in the cell using endoplasmic reticulum (ER) calcium resources. This process is modeled using modified Li-Rinzel model introduced by Tewari and Majumdar [22,43]: Here J chan denotes calcium influx from ER into cytosol through the IP3 receptor, J ERpump is the pumped calcium from intracellular space to ER and J ERleak indicates leaked calcium ions from ER into intracellular space (Details about f Slow can be found in the Appendix).
Presynaptic Ca 2+ ions can bind to Ca 2+ sensors of glutamate vesicles and contribute in vesicle release in the synaptic cleft [22,23,44]. The kinetic model for the calcium binding to the calcium sensor is governed by: Here, α is the rate of Ca 2+ association constant and β is the rate of Ca 2+ dissociation constant. γ, δ and β are constants which describe calcium-independent isomerization. X shows the Ca 2+ sensors on vesicles without any Ca 2+ bound. X(c i ) 1 denotes sensor with one Ca 2+ bound, X(c i ) 5 describes sensor with five Ca 2+ bound and finally Xðc i Þ In addition to evoked release due to the APs, the release of vesicles can occur spontaneously when there is no depolarization in the membrane. This type of release also depends on presynaptic Ca 2+ concentration and can occur when the Ca 2+ concentration is low in a non-depolarized membrane [44,45]. Spontaneous release rate due to Ca 2+ concentration is modeled by Nadkarni and Jung [45] using Poisson process and is expressed as follows: Where, a 1 = 50μM, a 2 = 5μM and a 3 = 0.85ms −1 are the Ca 2+ concentration at which λ is halved, slope factor of spontaneous release rate λ, and maximum spontaneous release rate, respectively. Values for these parameters are determined by Tewari and Majumdar such that, the frequency of spontaneous release is in the range of 1-3 Hz. Experimental findings indicate that this is a valid frequency range in the presence of astrocyte [22]. In general, glutamate release dynamics are depicted as follows [22]: Where, the fraction of ready to release vesicles is denoted by R, effective vesicles in synaptic cleft is represented by E, and inactive vesicles rate is denoted by I. τ inac = 3ms is the time constant of vesicle inactivation and τ rec = 800ms is the time constant of vesicle recovery. f r is the probability of vesicle release and has values (0, 0.5, 1). The total number of vesicles that can be released is assumed to be two vesicles. When two vesicles get ready to be released, the value of f r is 1, when one vesicle is ready to release value of f r is 0.5, and finally, when there are not any vesicles to be released f r is 0. The number of vesicles that get ready to be released can be calculated by stochastic simulation of Bollmann kinetic model for evoked release (by APs) [44] or by using Poisson function which refers to spontaneous release [22,23]. There is also an inactivation time of 6.34 ms for vesicle release (evoked or spontaneous release) [45]. Glutamate concentration in synaptic cleft can be altered due to some factors which are modeled by Tewari and Majumdar in the CA3-CA1 synapse. The fraction of effective vesicles, glutamate concentration in vesicles and docked vesicles and also degradation rate of glutamate in synaptic cleft can be considered as significant factors. Glutamate concentration can be represented by [22] Eq (15).
Here, g is the concentration of glutamate in the synaptic cleft, n v = 2 and g v = 60μM are the number of docked vesicles and glutamate concentration in vesicles respectively. g c = 10ms −1 is the degradation rate of glutamate which relates to neurotransmitter reuptake by astrocyte. Neurotransmitter reuptake is one of the major factors in the processing of addiction associated LTP formation and would be analyzed in this paper.

Astrocyte
Astrocytes release gliotransmitter, which can affect the synaptic transmission. Gliotransmitter release by astrocyte depends on the Ca 2+ concentration in astrocyte. The Released gliotransmitter can influence the presynaptic neuron through mGlurs, then postsynaptic neuron through NMDARs. In the proposed model, vesicle release dynamics by astrocyte are based on the Tewari and Majumdar's model [22]. Furthermore, It has been shown that during consumption of opioids, the astrocytic ability for reuptake is attenuated [28,29]. A decrease in astrocyte's glutamate reuptake ability is a primary cause for LTP induction that has been assessed during our simulation studies.
Released glutamate in the synapse, activates mGlurs of astrocyte, which leads to the augmentation of astrocytic Ca 2+ concentration. Ca 2+ Changes in astrocyte are because of IP 3 mechanisms. Ca 2+ oscillations in astrocyte is described by [46]: Where c Astro denotes Ca 2+ concentration in intracellular space of astrocyte, f Astro denotes the function which relates calcium currents to the calcium concentration. J chan,a is Ca 2+ flux from ER to the intracellular space due to IP 3 binding to ER's IP 3 R, J pump,a is the Ca 2+ removing from intracellular space by SERCA pump, J leak,a is the Ca 2+ leak from ER to intracellular space (detailed information about f Astro can be found in the Appendix).
Equations which describe vesicle release dynamics by astrocyte are [22] given by Eq (17).
þ Yðc Astro À c thresh Astro Þ:Pr a :R a ; Here R a denotes releasable vesicles in astrocyte, E a shows effective vesicles, I a is the inactive vesicles. Θ is Heaviside function, c thresh Astro ¼ 196:69nM is the threshold for calcium that can contribute to vesicle release, t a inact ¼ 3ms denotes to inactivation time constant and t a rec ¼ 800ms shows recovery time constant.
It is suggested that three independent gates are responsible for gliotransmitter release and each of the gates need a calcium ion to bind [22,47]. So, at least three calcium ions are needed for opening of the gates to permit vesicle release from astrocyte. This is described by [47] and depicted as Eq (18).
Here O j denotes the opening probability of the gate, k þ j is the opening rate of the gate, k À j is the closing rate of the gate and c Astro denotes calcium concentration O 3 denotes opening probability of a vesicle release site. Finally, Gliotransmitter dynamics in synapse is described by [22]: Here g a indicates gliotransmitter concentration, n v a ¼ 12 denotes the number of vesicles, g v a ¼ 20mM represents gliotransmitter concentration in one vesicle and g c a ¼ 10ðmsÞ À 1 shows clearance rate of glutamate.

Postsynaptic neuron
Spine of the postsynaptic neuron is a passive membrane which can produce EPSPs. Model of the spine consists of AMPA and NMDA receptors. Binding synaptic glutamate to AMPARs depolarizes the membrane potential. A depolarized membrane, potentially activates NMDARs in the presence of neurotransmitter and gliotransmitter. Activation of NMDARs allows the Ca 2+ ions to flow to the neuron. Ca 2+ entry precedes to CaMKII phosphorylation and NO synthesis which can retrograde to the presynaptic neuron. Phosphorylated CaMKII also lead to AMPAR insertion on the cell membrane, which is one of the signs of memory formation in synaptic scale. In the presence of morphine, NMDAR's current enhances because of morphine's modulation on channel conductance and morphine-induced inhibition of Mg 2+ blockage on NMDARs activity [14,31,35]. This enhancement of NMDAR current by morphine is dose dependent and was considered in the model. The CA1 pyramidal neuron's spine modeled using equations [48]: Here τ post denotes time constant of neuron membrane, V rest post shows neuron's membrane potential at rest, R m is the actual resistance of spin, I syn is synaptic current and can be described by: Where I AMPA and I NMDA are AMPA and NMDA receptors mediated currents respectively. AMPAR current considered using Destexhe's model [49]: Where V AMPA is reversal potential of the receptor, V post denotes membrane potential, m AMPA represents gating variable of AMPAR. g AMPA is the conductance of AMPAR and changes due to CaMKII phosphorylation process: Further information about AMPAR phosphorylation can be found in Castellani et all's papers [50,51].
Furthermore, in the proposed model, I NMDA modeled using Moradi et all's formulation with a slight modification: f NMDA represents a function that relates I NMDA to the receptor conductance (g NMDA ), gating variables of receptor (m NMDA ), Mg 2+ blocking (Mg) and postsynaptic membrane potential (V post ). Here in the proposed model, Mg 2+ blocking depends on morphine concentration (Morph) and membrane potential (V post ) through f Mg : Moreover, NMDAR gating variable is governed by: Here V NMDA is the reversal potential of NMDAR, α NMDA and β NMDA are opening and closing rate of receptor respectively, g NMDA is NMDA channels conductance which is described by: Which depends on the g VD (voltage dependent conductance) and morphine concentration through f g NMDA .
Enhancement of postsynaptic calcium is through AMPARs, NMDARs, and voltage-gated calcium channels and it descends is through calcium pumps. Postsynaptic calcium concentration can be described by [22]: Here, c post denotes postsynaptic calcium concentration which depends on AMPAR current, and is denoted by I AMPA , NMDA current is showed by I NMDA , voltage-gated calcium channels activation is described by I CaL and pumped calcium is defined by S pump .
Postsynaptic Ca 2+ concentration variation may lead to CaMKII phosphorylation which is governed by equations [52]: Finally enhancement of presynaptic calcium sensors sensitivity due to NOs is governed by a sigmoidal function: Here, P1 2 ¼ 40 mM is the half of total CamKII, k syt = 0.5% and k 1/2 = 0.4 μM are constants. Sensitivity of calcium sensors in Eq (12) enhances through this sigmoidal function.

Results
Simulation results of the proposed model are presented in this section. Simulations were implemented in Matlab 2014a software. Forward Euler method with fixed step size of 0.05 msec was used to solve the differential equations. Presynaptic neurons stimulated with a pulse train with 10 μA/cm 2 density (5 Hz with a duration of 4ms). Moreover, simulation time is 60 sec.

Mu-opioid receptor activation
In the first part, μORs activation was simulated. By applying morphine, μORs of interneurons is activated. Elevation in the morphine dose activates the receptor more. Bourinet et al. showed that by activating μORs, calcium influx through N-type calcium channels attenuated [20]. Inhibition of calcium channels depends on the amount of μORs activation [21]. Less calcium entry into the interneuron leads to less inhibitory neurotransmitter release. Therefore, the Inhibitory Postsynaptic Currents (IPSCs) in the CA3 pyramidal neuron decreases. Capogna et al. denoted that inhibition of IPSC by morphine in CA3 pyramidal neuron depends on morphine dose. This dependency can be modeled in a hill function form [19]. Simulation result for activation of μORs with the variety of morphine concentration has been shown in Fig 3. μORs activation leads to inhibition of interneuron calcium channels activation. So, the amplitude of IPSC in pyramidal neuron attenuates. This process is shown in Fig 4 in response to different doses of morphine. Fig 4A demonstrates that when there is not any morphine injection, all VGCCs are in an active state and the inhibitory current is at its maximum value. By applying 0.1 μM of morphine, VGCCs activation is attenuated to 75%, and IPSCs are decreased to 26 μA/cm^2 ( Fig  4B). A further increase in the dose of morphine to 1 μM precedes to 59% activation of VGCCs and IPSCs are attenuated to approximately 8 μA/cm^2 (Fig 4C).
Comparison of simulation results with the experimental findings for the activation of VGCCs and inhibition of the IPSC are denoted in Figs 5 and 6. During injection of 0.01 μM DAMGO (μOR agonist), calcium influx is decreased to 89% of its final value in the experiment [21]; simulation of the proposed model shows decrease in activated calcium channels to 95%. By increasing the dose of injected DAMGO to 0.1 μM, activated calcium channels attenuates to 87% in the experiment [21], and to 75% in the model. Further injection of DAMGO (1 μM) leads to more reduction in the calcium influx up to 66% in the experiment [21], and inhibition of activated channels decreases to 59% in the proposed model (Fig 5).
Inhibition of activated calcium channels in an inhibitory neuron can directly inhibit an inhibitory current on excitatory neurons based on Bertram's assumption. In this part, inhibitory current on the CA3 pyramidal neuron in simulation compared with experimental data in Fig 6. Different doses of morphine decrease the inhibitory current in the pyramidal neuron by 75%, 44% and 13% for 1, 0.1 and 0.01 μM injection of morphine in the experimental condition, respectively [19]. Simulation results for the inhibitory current show decrease by 92%, 57% and 5% With 1, 0.1 and 0.01 μM of morphine, respectively. One can see a similar pattern in simulation results and experiments. The model representing opioid receptor activation can relatively mimic the experimental results reported by Rhim et al. [21] and Capogna et al. [19].

Normal condition
To compare opioid-induced LTP with the typical situation, first, we simulate the model for the normal state, where there is no contribution of opioid receptors activation in the simulations. Stimulation of CA3 pyramidal neuron and interneuron with a pulse train leads to the generation of APs. APs activate N-type VGCCs and the Ca 2+ can enter into the neuron. Then, Ca 2+ binds to the presynaptic Ca 2+ sensors and triggers the glutamate release to the synaptic cleft ( Fig 7A). In addition, the effect of Interneuron activation is an augmentation of g GABA in the presynaptic pyramidal neuron which exerts a further inhibitory effect.
Synaptic glutamate can bind to astrocytic mGlurs and postsynaptic NMDA and AMPA receptors. By activation of astrocytic mGlurs, calcium concentration augments through IP3 production in astrocyte. Astrocytic Ca 2+ oscillations lead to the gliotransmitter release, which affects both presynaptic and postsynaptic pyramidal neurons. Gliotransmitter release has been shown in Fig 7B. Gliotransmitter affects presynaptic mGlurs and leads to IP3 production which is shown in Fig 7C. Production of presynaptic IP3 leads to the enhancement of calcium concentration in the synaptic bouton (Fig 7D). Synaptic glutamate can bind to postsynaptic AMPARs and NMDARs. Activation of AMPARs (Fig 7E) depolarize postsynaptic membrane potential and convey to remove of NMDARs Mg 2+ blocking, so NMDAR mediated current  (Fig 7F). Therefore, Ca 2+ influx to the postsynaptic neuron initiates. Postsynaptic Ca 2+ oscillations are shown in Fig 7G. Enhancement of Ca 2+ in the CA1 pyramidal neuron activates the CAMKII process which is able to augment AMPAR conductance and NO synthesis. In the Fig 7H, phosphorylated CaMKII is indicated by the red line and the threshold for the synthesis of NO with the green line. NO messenger can retrogradely diffuse from postsynaptic neuron to the presynaptic neuron and increases the sensitivity of Ca 2+ sensors in the CA3 pyramidal neuron. Subsequently, there would be a further glutamate release in the synaptic cleft. This process acts as a positive feedback to inducing and maintaining the LTP in the synapse. Overall, this simulation shows a typical learning process in synaptic level.

Pathological condition
Injection of 1 μM morphine activates presynaptic and postsynaptic μORs. Activation of presynaptic μORs inhibits the Ca 2+ entrance into the interneuron by 40%, which is described in Fig 5. Therefore, IPSC of CA3 pyramidal neuron decreases to 8 μA/cm^2 due to fewer GABA releases by interneuron. Reduction of IPSC leads to more neuronal excitability. Consequently, synaptic glutamate increases and stimulates the astrocyte and postsynaptic neuron more than its normal condition. Synaptic glutamate has been shown in Fig 8A. As one can see, Computational model opioid synpatic plasticity amplitude, frequency and the total amount of released glutamate by the presynaptic neuron is significantly increased. This increase is about six times larger than its normal condition and may lead to more excitation of both postsynaptic neuron and astrocyte. Furthermore, gliotransmitter release by astrocyte is enhanced by 250% of its normal state which is described in Fig 8B. Gliotransmitter increase is due to further activation of astrocytic mGlurs and enhancing calcium oscillations of the astrocyte. Binding gliotransmitter to presynaptic mGlurs acts as a positive feedback and results in IP3 production and increase of presynaptic calcium, which are denoted in Fig 8C and 8D, respectively. As it has been shown in the Fig 8C, IP3 production is enhanced in the presynaptic neuron, which can suggest the increase of calcium in synaptic bouton. Furthermore, since glutamate reuptake is attenuated in the presence of opioids, we modeled it by decreasing the reuptake ability of glutamate by astrocyte down to 50% of a typical situation based on experimental observations. On the other hand, the presence of opioids affects NMDAR activation by inhibiting Mg 2+ blocking and by increasing conductance of NMDA receptors. Due to experimental research, Mg 2+ blocking of NMDARs modified in which δ (the relative electrical distance of the binding site of Mg 2+ from the outside of membrane) increased by 11% and K 0 (dissociation constant at 0 mv) augmented by 4.8 fold of its normal value. Besides that, morphine can increase NMDARs voltage independent conductance by 15% averagely, by Protein Kinase C (PKC) mechanisms, which has been applied in the model.
Generally, in the pathological conditions, postsynaptic calcium's entry to the neuron has been denoted in Fig 8G, which shows a significant increase in comparison to the normal state. This amplification of calcium oscillations is due to the enhancement of NMDARs activation (Fig 8F). NMDAR related current increases because it receives more glutamate from the presynaptic neuron and more gliotransmitter from astrocyte. Besides that, activation of the μORs enhance conductance of the NMDARs through the previously mentioned mechanisms. Calcium entry to the postsynaptic neuron causes phosphorylation of CaMKII which is shown in Fig 8H. Phosphorylation of CaMKII is followed by NO synthesis. Synthesized NO retrogrades to the presynaptic neuron and enhances the binding ability of glutamate vesicles to calcium Computational model opioid synpatic plasticity ions. Besides that, enhancing calcium entry to the postsynaptic neuron also leads to an increase in AMPAR conductance by approximately 30% in pathological condition. Generally, with the assumption that enhancement of the CaMKII phosphorylation, AMPARs conductivity, glutamate, and gliotransmitter are the main traces of the induction of LTP in the synapse, morphine injections appear to lead to pathological memory.
It has been shown that opioid presence leads to LTP production. Fig 9 compares normal and pathological conditions in two block diagrams. Under normal conditions, the interneuron regulates the activity of the presynaptic pyramidal neuron, and inhibitory current is significant. Also, astrocytic transporters work normally (Fig 9A). In the pathological conditions, the inhibitory current decreases and the amount of neurotransmitters and gliotransmitters increases. Also, by applying morphine, the sensitivity of the postsynaptic NMDARs to the Computational model opioid synpatic plasticity transmitters enhances. Therefore, more calcium enters into the neuron. The increase in calcium conveys to CaMKII phosphorylation and thus NO is synthesized. Also, AMPAR's conductance increases due to CaMKII phosphorylation process (Fig 9B).

Comparing the effect of different factors on pathological LTP induction
In the previous simulations, it has been demonstrated that morphine injection conveys to pathological LTP induction. In this section, the effect of the different factors on the induction of pathological LTP is considered and significant parameters are specified. Essential elements can be seen as drug targets to prevent the formation of pathological memory in the use of opioids. The presence of 1 μM morphine in synapse can enhance CaMKII phosphorylation process, synaptic glutamate, gliotransmitter release, and AMPAR conductance. These variables are considered as monitoring factors for morphine effects on the synapse. Monitoring elements can be seen as a trace of LTP induction. Fig 10 shows  In Simulation 1 morphine affects presynaptic, postsynaptic and astrocyte. Therefore, presynaptic and postsynaptic μORs are in the active state, and astrocytic glutamate transporters coefficient is half of the normal conditions, which indicates the effect of morphine on astrocyte.
Comparing these results with normal conditions (Simulation 13) infers that the average amount of synaptic glutamate is five times higher than normal. The average amount of gliotransmitter is 4.5 times the normal state. The mean amount of phosphorylated CaMKII and AMPARs are 4 and 1.3 times higher than the normal level. Overall, in the presence of morphine, all of the monitoring factors are high, which can be indicative of LTP induction. Redline shows the median of the signals, '+' sign represents the mean value of each signal, blue line demonstrates that 25%-75% of data are included in this range, and finally black line shows the range that 9%-91% of the dada are included. Panel A shows the synaptic glutamate, Panel B denotes the gliotransmitter, Panel C describes the phosphorylated CaMKII, and Panel D represents the Phosphorylated AMPARs. In the first simulation, all of the synaptic components are affected by morphine and the simulation of 13 th indicates normal mode. In the second simulation, the effect of morphine on presynaptic neurons has been neglected. In the third simulation, postsynaptic μORs are assumed to be inactive. In the fourth simulation, morphine only affects astrocyte. In the fifth simulation, the release of the gliotransmitter is considered zero. In the sixth simulation, it is assumed that despite the effect of morphine on all the components of the synapse, astrocytic GLTs act normal. In the seventh simulation, under conditions where morphine affects all of the synaptic components, the activity of astrocytic transporters is stimulated to be 50% greater than normal. In the eighth simulation, the activity of astrocytic mGlurs is considered zero. The ninth simulation shows that astrocytic mGlurs and postsynaptic μORs are not active. In the tenth simulation, astrocytic transporters are stimulated 50% more than normal, and astrocytic mGlurs have been deactivated. In the eleventh simulation, by stimulating astrocytic transporters, NMDA receptors have a normal synaptic effect. In the 12 th simulation, astrocytic mGlurs have been deactivated, astrocytic transporters are stimulated, and NMDARs act as normal. https://doi.org/10.1371/journal.pone.0193410.g010 Computational model opioid synpatic plasticity Simulation 2. In this simulation, we want to show how the disinhibitory effect of morphine affects monitoring factors. Here, we disable the presynaptic opioid receptors and morphine only affects the postsynaptic μORs and astrocytic GLT1s.
Simulation results denote that despite a significant reduction in synaptic glutamate, other monitoring factors are still high and pathological LTP is still present. This means that, despite the importance of disinhibitory mechanism in inducing LTP, preventing it cannot block pathological memory.
Simulation 3. In this simulation, the postsynaptic μORs activity is inhibited. however, presynaptic neurons and astrocyte are under the influence of morphine. The results indicate that, despite the attenuation of CaMKII phosphorylation and the reduction of AMPAR conductance, the amount of glutamate and gliotransmitter are still highand the probability of formation of pathologic memory is high.
Simulation 4. In this section, morphine is thought to only affect astrocyte. Therefore, presynaptic and postsynaptic μORs are inhibited. In this case, the synaptic glutamate is approaching a level that is normal. In addition, phosphorylation of CaMKII and AMPARs has been significantly reduced. But it seems that due to the significant release of glutamate in the synapse. LTP occurs. This simulation describes that astrocyte plays a fundamental role in the induction of pathological LTP, so it is necessary to carefully examine its role in the following sections.
Simulation 5. In this simulation, it has been assumed that morphine affects all the components of the synapse, but astrocyte cannot release the gliotransmitter. Here astrocyte does not have any effect on presynaptic and postsynaptic neurons. The results suggest that inhibiting gliotransmitter release cannot disrupt the pathological memory formation. Simulation 6. In this simulation, the role of astrocytic transporters in pathological memory induction has been investigated. Therefore, under the conditions of the first simulation, it is assumed that astrocyte transports are active as normal and their activity is not affected by the presence of morphine. The results show that normal activity of transporters cannot prevent the formation of pathologic memory.
Simulation 7. In this section, when morphine affects all the components of the synapse, the activity of astrocytic transporters enhanced by 50% more than their standard value. In this case, glutamate and gliotransmitter values have gone to normal, and CaMKII phosphorylation has been significantly reduced. This simulation shows that up-regulation of astrocytic transporters can be used as an agent to prevent the pathological LTP formation.
Simulation 8. In this simulation, the effect of astrocytic mGlur has been investigated. To this purpose, astrocytic mGlurs are inhibited in a situation where morphine affects all of the synaptic components. It can be seen that in this case, with the exception of the gliotransmitter which is in the normal range, the rest of the monitoring factors are high.
Simulation 9. In this simulation, postsynaptic opioid receptors and astrocytic mGlur are supposed to be inhibited. The results state that it is still possible to develop pathological memory.
Simulation 10. In this simulation, astrocytic mGlurs is inhibited and its transporters are stimulated 50% more than normal. Opioid receptors for presynaptic and postsynaptic neurons are still active. Results indicate that in spite of a significant drop in glutamate and gliotransmitter, the rest of the monitoring factors were not significantly reduced.
Simulation 11. In this simulation, postsynaptic opioid receptors have been disabled and astrocytic GLTs are stimulated 50% more than normal. The results show that in the presence of morphine, all of the monitoring factors are in normal range.
Simulation 12. In this simulation, postsynaptic opioid receptors and astrocytic mGlurs are supposed to be inhibited and astrocytic transporters are stimulated 50% higher than normal. In this case, monitoring factors are found to be close to the normal value, and even take less than normal values. Simulation 13. This simulates the normal state in which it is assumed that morphine is not present in the synapse.

Discussion
This paper presented a computational model in CA3-CA1 region of the hippocampus at the synaptic level. Here we extended and improved previously introduced models by different researchers [22,23,34,41,42,45,46,[48][49][50][51][52][53] to represent opioid-induced synaptic plasticity in the hippocampus. To best of our knowledge, this is the first model that has the ability to be used for computational analysis of opioid-induced synaptic plasticity. Using computational models allows researchers to analyze the effect of various parameters simultaneously without concerning about the interfering effects of different medical agents in experiments.
Our model adequately shows that opioids presence can induce LTP in affected synapses. Here, it is assumed that amounts of synaptic glutamate, gliotransmitter, phosphorylated CaM-KII and AMPARs are footprints of LTP induction in the synapse. In fact, these factors are considered to be parameters which can lead to LTP induction, and we called them as monitoring factors. To discuss the results of the simulations, the normalized values for monitoring factors are depicted in Fig 11. As previously stated, the first simulation shows that morphine affects all the components of the synapse and the thirteenth simulation indicates the normal mode. It has been shown in experiments that the presence of morphine increases the density of AMPARs, synaptic glutamate and calcium signaling in affected neurons [9,54]. These changes, which describe the main features of pathological LTP and memory formation, are seen in the results of the first simulation.
In the second simulation, presynaptic opioid receptors are inhibited. As a result, release of glutamate is normal, but the traces of pathological LTP is still existed. The results of this simulation indicate that increasing the activity of postsynaptic NMDARs and the effect of morphine on astrocyte are factors that contribute to the formation of pathologic memory.
In the 3 rd simulation, postsynaptic μORs have been inhibited, indicating that postsynaptic NMDARs exhibit normal activity. In this case, the monitoring factors have significant amounts. It infers that injection of NMDAR antagonist or inhibition of postsynaptic opioid receptors is effective in preventing the induction of pathologic LTP, but it is not enough.
In fact, experimental results have also shown that NMDARs play an essential role in inducing opioid-related neural plasticity [38,[55][56][57]. Furthermore, it has been shown that NMDAR antagonist injections are effective in relapse weakening; however, it cannot completely eliminate drug-seeking behavior [58,59]. Based on the results we have obtained in the simulations, it is true that NMDARs are the main actors in the formation of pathological memories, but the role of other parameters should not be ignored.
In the fourth simulation, only astrocyte is affected by morphine, and presynaptic and postsynaptic opioid receptors are inhibited. In this case, glutamate is found to be in normal range, but the amounts of CaMKII and AMPAR phosphorylation are still significant. As a result, it can be concluded that, regardless of the role of astrocyte in synapse, the pathologic LTP cannot be controlled. Some researchers believe that astrocyte dysfunction can lead to addiction [28,30]. Thus, it can be concluded that without the manipulation of astrocyte function, we may not be able to inhibit the pathological LTP induction.
In addition, in the 5 th and 6 th simulations, it was observed that preventing the release of gliotransmitter and having normal GLTs cannot influence the formation of pathological LTP.
In the 7 th simulation, astrocytic transporters have been stimulated by 50% more than normal. It can be seen that despite the attenuation of synaptic glutamate and gliotransmitter, other monitoring factors are still high. So, stimulating astrocytic GLTs can be effective, but is not sufficient for preventing pathological LTP. Now the important question is: is it possible to up-regulate astrocytic glutamate transporters empirically? The answer is yes. It has recently been shown that up-regulation of astrocytic transporters is possible and can even undermine the induction of LTP in the synapse [60].
In the 8 th simulation, astrocytic mGlurs are inhibited. This mode describes that the agonistdependent IP3 production is inhabited in astrocyte. According to the results obtained, this mode alone cannot be used as an effective parameter.
In the 9 th simulation, astrocytic mGlurs and postsynaptic μORs are simultaneously inhibited. In this case, synaptic glutamate is high, phosphorylated CaMKII is almost halved, and the amount of phosphorylation of AMPARs and released gliotransmitter are in normal range. It can be noticed that blocking the agonist-dependent IP3 production and the use of NMDAR antagonist can undermine the pathological LTP. In this case, the only remaining problem is to reduce the amount of synaptic glutamate that can also regulate CaMKII phosphorylation.
In the 10 th simulation, activity of astrocytic mGlurs is inhibited, and GLTs are stimulated to be 50% more than normal level. In this case, glutamate and gliotransmitter are in the normal range, and the other monitoring factors are almost halved. Therefore, it can be concluded that manipulation in the activity of astrocyte can be effective in preventing the formation of pathological LTP.
The 11 th simulation describes the case in which astrocyte transporters are stimulated, and postsynaptic μORs are inhibited. In this case, all monitoring factors fall into the normal range. It means that stimulating GLTs and injection of NMDAR antagonists can stop the formation of pathological LTP.
In the 12 th simulation, in addition to inhibiting the activity of postsynaptic opioid receptors and stimulating GLTs, we also disable astrocytic mGlurs. In this case, monitoring factors are alose to normal and even some of them have values that are less than normal. It can be concluded that in this case the possibility of the formation of pathological memory is completely eliminated.
Overall, by summerizing all the simulations, it can be suggested that an optimal method for coping with the formation of pathological memory can be the stimulation of astrocyte GLTs and the injection of NMDAR antagonists.
In this paper, some assumptions are made to simplify the modeling process. Dynamics has not been considered for the release of inhibitor neurotransmitter from the interneuron to the pyramidal neuron in the presynaptic layer. Such an approach has been implemented in some computational models. It is also assumed that interneuron and astrocyte do not interact. Investigation of this interaction and its role in the induction of the synaptic plasticity can be considered in the future works. Also, in this paper, as with some computational models, the diffusion process for glutamate and gliotransmitter is not considered. It seems that focusing on this process in models that deal with a large number of neurons (network level) can be effective, but in our model, it will not have much impact on the results. On the other hand, we assume that some factors are footprints of LTP induction in the synapse, such as the concentration of glutamate and gliotransmitter, and phosphorylation of CaMKII and AMPARs. Considering the effect of LTD in synapse can extend the model. This goal may be accomplished through the use of the STDP process as a learning rule in synapse.

Conclusion
Opioids such as morphine can lead to the induction of LTP at the synaptic level, which may result in an abnormal memory formation process leading to addiction. This pathological memory can be seen as the key trigger for relapse in addicted patients. Suppression of LTP induced by opioids may lead to diminishing relapse in withdrawal period in addicted subjects. Since changing different parameters and evaluating various types of effects in memory formation process at the synaptic level can be difficult and even practically impossible in experimental conditions, a detailed mathematical model can help us to overcome empirical obstacles and difficulties to some extent. This paper, presented a biophysically suggested and supported computational model for opioid-induced memory formation process. The proposed model is at the synaptic level and the role of astrocyte is considered as the third component of synapse. The proposed model is capable of describing the pathologic LTP induction process and can be helpful in addiction-related research. Here, we found that astrocytic transporters play a fundamental role in the induction of pathological drug-related memories. Simulation results show that up-regulating astrocytic transporters, omitting astrocytic mGlurs and applying NMDARs antagonist can completely prevent pathological LTP. We also suggest that the up-regulation of GLTs and injection of NMDAR antagonist can be used as an optimal way to prevent the occurrence of pathological LTP. Overall, we suggest that astrocyte has a prominent role in addiction-related memories and its function must be considered precisely for preventing abnormal memories. This approach may help researchers to inhibit relapse in the future. astrocyte can be written as follows: Hillðc 2 a ; K PLCd Þ þv 3K Hillðc 4 a ; K D ÞHillðp a ; K 3 Þ À r 5p a p a dh a dt ¼ a h a ð1 À h a Þ À b h a h a ðA-9Þ The first term of above equation denotes agonist-dependent (glutamate) IP 3 production; the second term shows agonist-independent IP 3 production which is due to PLCδ signaling and modulated by Ca 2+ , the third term shows IP 3 degradation rate by IP 3 − 3k which is Ca 2+ dependent and the last term incorporates degradation of IP 3 by IP − 5P. a h a and b h a denote opening and closing rate of h a respectively. Other parameters in mentioned equation can be determined using equations: m 1;a ¼ Hillðp a ; d 1 Þ; n 1;a ¼ Hillðc a ; d 5 Þ; Hillðx n ; KÞ ¼ x n x n þ K n ; a h a ¼ a 2 d 2 Further information about the parameters and their values can be found in De pitta's paper [46].
Postsynaptic neuron. AMPAR current governs by the following equations: I AMPA ¼ f AMPA ðg AMPA ; m AMPA ; V post Þ ¼ g AMPA m AMPA ðV post À V AMPA Þ ðA-11Þ Where V AMPA is reversal potential of the receptor, V post denotes membrane potential, m AMPA represents gating variable of AMPAR. Since we assume that both presynaptic and astrocytic glutamate have effects on postsynaptic neuron, so gating variable of AMPAR can be described by: dm AMPA dt ¼ a AMPA ða 1 g pre þ a 2 g ast Þð1 À m AMPA Þ À b AMPA m AMPA ðA-12Þ Here α AMPA and β AMPA are opening and closing rate of the receptor, g pre and g ast denote glutamate concentration released by presynaptic neuron and astrocyte respectively, a 1 = 0.42 and a 2 = 0.01 are constants.
In Eq (A-11) AMPAR channel's conductance (g AMPA ) can be varied due to postsynaptic CaMKII. This process modeled using modified equations [50,51] Here, the activity of protein kinase and protein phosphatase assumed to be a hill function of CaMKII and have been shown by EK and EP respectively [51].

ðA-14Þ
Here, A, A p 1 , A P 2 and A P 2 P 1 are reaction substrates and products. The conductance of AMPAR can be determined using the below equation:
Furthermore, in the proposed model, I NMDA modeled using Moradi et all's formulation with a slight modification: Here, Mg denotes Mg 2+ blocking and is governed by: Moreover, NMDAR gating variable is governed by: dm NMDA dt ¼ f m NMDA ða NMDA ; b NMDA ; g pre ; g Astro Þ ¼ a NMDA ðk 1 g pre þ k 2 g Astro Þð1 À m NMDA Þ À b NMDA m NMDA ðA-18Þ Here V NMDA is the reversal potential of NMDAR, α NMDA and β NMDA are opening and closing rate of receptor respectively, g pre and g Astro are glutamate concentrations released by presynaptic pyramidal neuron and astrocyte respectively. g NMDA is NMDA channels conductance which is described by: Here g VI is the channel conductance independent of potential and g VD is the channel conductance dependent to the potential which can be written by: dg VD dt ¼ ðg VD;1 À g VD Þ=t g g VD;1 ¼ kðV post À V 0 Þ

ðA-20Þ
Where g VD,1 is the final value of g VD with time constant of τ g and has a linear relation with membrane potential with constant of k. Other parameters and values listed in S1 Table. catalytic constants respectively. Then, phosphorylated CaMKII can be described by: