1. Introduction
Opioid consumption like other drugs can lead to behavioral changes (Contet, Kieffer, & Befort, 2004; Williams, Christie, & Manzoni, 2001). Some researchers believe that these changes are due to pathological memory formation which can be called memory of addiction (Dong & Nestler, 2014; Kelley, 2004; Nestler, 2001). Hippocampus seems to be one of the targets subject to the plastic changes due to drug abuse (Bao et al., 2007; Borjkhani, Mahdavi, & Bahrami, 2014; Capogna, Gähwiler, & Thompson, 1993; Caudle & Chavkin, 1990; Ga, 1980; Heidari et al., 2013; Hosseinmardi, Azimi, Fathollahi, Javan, & Naghdi, 2011; Hosseinmardi, Fathollahi, Naghdi, & Javan, 2009). 
At the synaptic level, drug-induced changes might be due to Long-Term Potentiation (LTP)/ Long-Term Depression (LTD) induction. Therefore, analyzing LTP/LTD induction in synapses under the presence of opioids can help researchers to unveil the formation mechanisms of this pathological memory by drugs (Kauer & Malenka, 2007; Lüscher & Malenka, 2011). 
Opioids affect hippocampal neurons mainly by enhancing glutamate release in presynaptic neurons (Akaishi, Saito, Ito, Ishige, & Ikegaya, 2000; Williams et al., 2001) and increasing postsynaptic receptors activities (Kauer & Malenka, 2007; van Huijstee & Mansvelder, 2015; Williams et al., 2001). All the modifications caused by opioids result in LTP induction in affected synapses (Bao et al., 2007; Wolf, 2003). It has been shown that one of the causes of LTP induction is Calmodulin-dependent protein Kinase II (CaMKII) phosphorylation process (Lisman, Yasuda, & Raghavachari, 2012; Sanhueza & Lisman, 2013). 
On the other hand, Ca2+ concentration is the main mediator in this process. Elevation of Ca2+ concentration in postsynaptic neuron triggers CaMKII phosphorylation process. Ca2+ can enter into the neuron through Voltage-Gated Calcium Channels (VGCCs), N-Methyl-D-Aspartate Receptor (NMDAR) and Ca2+ permeable α-Amino-3-hydroxy-5-methyl-4-isoxazole Propionic Acid Receptors (AMPARs). Thus, these channels and receptors contribute in CaMKII phosphorylation process and LTP induction. 
It has been shown that blocking hippocampal VGCCs can prevent opioids dependency (Mishra, Barik, & Ray, 2017). Furthermore, opioid effects on NMDARs and AMPARs were verified in different research studies (Capogna et al., 1993; Chen & Huang, 1992; Chen & Marine, 1991; Ga, 1980; Kauer & Malenka, 2007; Nestler, 2013; Peters & De Vries, 2012; Williams et al., 2001). Most of the research and findings in this field are based on experimental observations. 
However, computational models offer new tools and approaches for investigating neurophysiological basis of addiction memory formation process. In general, computational modeling of biological systems helps us to study consistency of different experimental findings and to predict the behavior of a given biological system under different constraints. That is why computational modeling has recently attracted the interests of many researchers in the field of neuroscience (Amiri, Bahrami, & Janahmadi, 2012; Amiri, Montaseri, & Bahrami, 2011; Li et al., 2016; Tewari & Majumdar, 2012b; Volman, Bazhenov, & Sejnowski, 2012). 
Therefore, in this paper, using experimental findings we introduce a minimal computational model to assess opioid-induced LTP in CA1 region of the hippocampus. This model can also be used for analyzing the effects of other drugs by minor modifications. In the first part of the paper, opioids effect on hippocampal neurons will be discussed. The second part is dedicated to introduce mathematical modeling approach. The third part presents simulation results. Discussion and conclusion will be presented in the final part. 

2. Methods 
2.1 Neurophysiological background
According to many research studies, opioid-induced modifications in glutamatergic synapses are responsible for pathological memory formation related to the addiction (Harris, Wimmer, Byrne, & Aston-Jones, 2004; Herman et al., 2003; Mameli, Bellone, Brown, & Lüscher, 2011; van Huijstee & Mansvelder, 2015). Since Long-Term Potentiation (LTP)/ Long-Term Depression (LTD) induction can be considered as signs of memory formation process at the synaptic level, one can also suggest that addiction memory formation might occur through the induction of LTP/LTD. 
It has been shown that opioids induce LTP in hippocampal synapses by two primary mechanisms: 1. Enhancing glutamate release through disinhibitory mechanism (Capogna et al., 1993; Caudle & Chavkin, 1990; Cohen, Doze, & Madison, 1992; Rezai, Kieffer, Roux, & Massotte, 2013); and 2. By manipulating postsynaptic NMDARs current (Chen & Huang, 1992; Chen & Marine, 1991; Kow, Commons, Ogawa, & Pfaff, 2002; van Huijstee & Mansvelder, 2015). Opioids modify NMDARs current by increasing conduction of NMDARs and inhibiting Mg2+ from blocking NMDA receptors (Chen & Huang, 1992; Chen & Marine, 1991; Herman et al., 2003). 
Consequently, these changes in synaptic transmission result in more calcium flux into the postsynaptic neuron (Kow et al., 2002; Przewlocki et al., 1999). Enhancement of postsynaptic calcium can activate CAMKII mechanisms (Herman et al., 2003; Trujillo, 2002). Phosphorylation of CaMKII leads to LTP induction and can be seen as molecular memory formation (Kauer & Malenka, 2007; Lisman et al., 2012; Sanhueza & Lisman, 2013). Furthermore, enhancement of postsynaptic calcium can lead to phosphorylation of AMPARs (Castellani, Quinlan, Bersani, Cooper, & Shouval, 2005; Castellani, Quinlan, Cooper, & Shouval, 2001; Kauer & Malenka, 2007). 
Also, it has been shown that VGCCs are involved in opioids dependency (Bongianni, Carla, Moroni, & Pellegrini‐Giampietro, 1986; Mishra et al., 2017). To investigate the mechanism of a false LTP induced by opioids, we developed a computational model of a postsynaptic neuron in CA1 region. In this model, we considered all important aspects mentioned from these neurophysiological studies. This model will be described in the next part (Figure 1). 

2.2. Computational modeling approach
To analysis opioid-induced synaptic plasticity in CA1 region of the hippocampus, the spine of a CA1 pyramidal neuron is modeled by Equation (1): 

 

This Equation denotes a minimal model to produce Excitatory Postsynaptic Potentials (EPSPs) based on Tsodyks and Markram’s passive membrane model (Tsodyks & Markram, 1997), where τpost shows time constant of neuron membrane, Vrestpost denotes neuron’s membrane potential at rest, and Rm denotes actual resistance of spin. The synaptic current denoted by Isyn can be described by:

 

Here, IAMPA and INMDA are AMPAR and NMDAR currents, respectively. Model of AMPAR current is based on Destexhe’s model (Destexhe, Mainen, & Sejnowski, 1998):

 

Where VAMPA refers to reversal potential of the receptor, Vpost is the membrane potential, and mAMPA denotes gating variable of AMPAR. Gating variable of AMPAR can be represented by the following Equation (Destexhe, Mainen, & Sejnowski, 1998):

 

Here opening and closing rate of the receptor are shown by αampa and βampa, respectively, gpre shows glutamate concentration, gAMPA described in Equation (3) denotes AMPAR’s channel conductance that can be enhanced due to CaMKII phosphorylation process by using the following Equation:

 

Where gAMPA0=0.4 nS is the initial value for AMPAR’s conductance, Phalf=40 µM and khalf=0.4 µM are constants which enables AMPAR conductance to increase through CaMKII phosphorylation process. Model of the NMDAR current can described by the following Equation (Destexhe, Mainen, & Sejnowski, 1998; Moradi et al., 2013):

 

Where VNMDA shows the reversal potential of NMDAR, opening and closing rate of receptor are shown by αNMDAand βNMDA respectively, Op refers to opioids concentration in the synapse, and gpre shows glutamate concentration. NMDA channels conductance (gNMDA) is represented by the following Equation:

 

, where gVI and gVD are voltage-independent and voltage-dependent conductance. The latter conductance is described by: 

 

Here, gVD,∞ is the final value of gVD, τg denotes time constant, gVD,∞ , and Vpost has a linear relation with constant of K. Other parameters and their values have been listed in the Appendix.
Calcium can enter into the neuron through AMPARs, NMDARs, and VGCCs. Also, calcium pumps are responsible for transferring the calcium out of the neuron. So, calcium concentration inside the neuron can be represented by the following Equation (Tewari & Majumdar, 2012b: 

 

Here, Cpost shows postsynaptic calcium concentration, IR denotes voltage-gated calcium channels activation, and Spump refers to pumped calcium. Here, η=0.012 and y=0.06 demonstrate the amount of calcium that can enter into the neuron through AMPARs and NMDARs. Zca=2 denotes calcium valence, F=96487C/mo l is the faraday’s constant, and Vspine=0.9048 µM3 is the volume for dendrite spin. Furthermore, bt=200 µM shows total endogenous buffer concentration and Kendo=10 µM denotes endogenous buffer calcium affinity. Pumped calcium is represented by the following Equation (Tewari & Majumdar, 2012b): 

 

Here, Ks=100/s is the maximum efflux rate of calcium pump, Crestpost=100nM is the rest value for postsynaptic calcium concentration, and VGCCs activation is represented by the following Equation (Tewari & Majumdar, 2012b):

 

, where gR=15ps is the conductance of calcium channel, B(N, Popen) denotes a random variable with a binomial distribution which shows the number of opened channels, and VR=27.4mv denotes the reversal potential. 
Phosphorylation of CaMKII due to calcium concentration can be described by the following Equation (Zhabotinsky, 2000):

 

 

, where Pi shows the concentration of the i-fold phosphorylated CaMKII, ep refers to the PP1 concentration which is not bounded to l1P and may demonstrate active protein phosphatase, the total concentration of PP1 is shown by ep0=0.1 µM, free l1P is denoted by I, and free l1 concentration is shown by I0=0.1 µM. Also K3=1/ µMs and K3=10-3/s are association and dissociation rate constant of PP1-l1P complex, respectively. vcaN=2/s refers to the rate of l1P dephosphorylation due to calcineurin (CaN), VPKA=0.45 µM/sis the phosphorylation rate of l1 due to the protein kinase A (PKA), and KH2=0.7 µM is the calcium activation Hill constant of CaN. 
Phosphorylation (V1), auto-phosphorylation (V2), and dephosphorylation (V3) rates can be represented by the following equations:

 

, where k1=0.5/s is the l1 dependent regulation rate of PP1 and KH1=4 µM is the Hill constant of CaMKII for calcium activation. KM=20 µM and are K2=10/s the Michaelis and catalytic constants, respectively. Finally, phosphorylated CaMKII can be presented by:

 

In the next part, simulation results are presented. 

3. Results
The mathematical modeling approach simulation results are presented in this section. Simulations were implemented in Matlab 2014a. To solve differential equations, the forward Euler method with constant step size of 0.05 ms is used. Synaptic glutamate is assumed as a periodic signal with a frequency of 5 Hz which each pulse has a duration of 4 ms. The amplitude of this signal is 0.2 mM and the simulation time is 10 s. In all Figures, 1 second of the simulation is depicted to clarify dynamics of the signals. 

3.1. Normal and pathological conditions
One of the features of opioid-induced synaptic plasticity is the induction of LTP where normal LTP cannot occur (van Huijstee & Mansvelder, 2015). Therefore, stimulation signal (which is synaptic glutamate) is considered as a pulse train with a frequency of 5 Hz and amplitude of 0.2 mM. Choosing these values are based on experimental observations used by Tewari and Majumdar’s computational models (Tewari & Majumdar, 2012b; Tewari & Majumdar, 2012a). 
 

 

As it is shown, synaptic glutamate activates AMPARs and NMDARs. By activation of these receptors, calcium enters into the neuron through VGCCs, NMDARs, and AMPARs. Activation of calcium pumps decreases calcium concentration in the neuron. Elevation of calcium may lead to activation of CaMKII mechanism which enhances AMPAR conductance. 

 

 

 

The plain black bar is the phosphorylated CaMKII, and the green bar with black patterns is the phosphorylated AMPAR. Horizontal axes with “normal” label shows simulation in normal condition, “Opioid” label shows opioid injected condition, “CK=0” stands for the simulation in which opioid effect on NMDAR's conductance is assumed as 0, “CMg=0” describes the state that opioid effect on Mg2+ blocking is assumed as 0 and “Icr=0” refers to the case that current mediated by VGCC is assumed as 0.
Applying synaptic glutamate (Figure 2a) leads to AMPARs activation (Figure 2b) which depolarize membrane potential. Depolarized membrane allows NMDARs activation by removing Mg2+ blocking (Figure 2c). Since high-voltage calcium channels are considered in the model, calcium channels activation is 0 which is described in Figure 2d. Ions of Ca2+ can enter into the neuron through AMPARs and NMDARs activation. So, Ca2 concentration enhances in the neuron which is shown in Figure 2e. Entered Ca2+ can engage in CaMKII phosphorylation process, however, due to lack of enough calcium, this process is not observed here (Figure 2f). Based on simulation results, none of the AMPARs can phosphorylate (Figure 2g). Generally, this simulation can describe the normal condition which there is not any opioid’s effect in the synapse. In this simulation, LTP does not occur in the synapse due to lack of phosphorylated CaMKII and AMPARs. 
Applying 1 μM of opioid increases NMDARs activation through inhibition of Mg2+ blocking and enhancement of receptors conductance. NMDAR mediated current enhances in this situation (Figure 3c). Thus, more calcium enters into the neuron (Figure 3e) which triggers CaMKII phosphorylation process (Figure 3f). Furthermore, more activation of NMDARs results in VGCCs activation. Figure 3d shows calcium current mediated by VGCCs. Exceeding phosphorylated CaMKII from a predetermined threshold can phosphorylate AMPARs. In fact, AMPAR conductance increases which is shown in Figure 3g. Thus more current is mediated by AMPARs (Figure 3b) that can work as a positive feedback for enhancement of calcium concentration in the neuron. 

3.2. Comparing different factors on pathological LTP induction
Simulation results show that injection of opioid can induce LTP when normal LTP does not occur. Now we want to assess which components of the model are responsible in inducing such a pathological LTP. Based on these simulations we can suggest a way to prevent this type of pathological LTP. 

 

 

 

Further simulations were performed to analyze different parameters involved in inducing pathological LTP. Normalized value for phosphorylated CaMKII and AMPARs are calculated in each simulation. The values are calculated by using the following Equation:

 

Where  YNormalized refers to the normalized value, Y is the parameter value, also YMax and YMin denote maximum and minimum values of the parameter, respectively. This normalization approach results in normalized values within [0-1] interval. 
Opioid-induced LTP were analyzed with 3 different assumptions: 1. Opioid does not modify NMDARs conductance; 2. Opioid does not have any infulence on Mg2+ blocking; and 3. VGCCs mediated current is 0. Simulation results according to the assumptions are shown in Figure 4. In normal condition, phosphorylated CaMKII and AMPAR is 0 which specifies there is not any step foot of LTP in the synapse. Presence of opioid leads to phosphorylation of CaMKII and AMPAR, and normalized values are in the high state. This simulation demonstrates the presence of LTP in the synapse. 
It can be seen that LTP does not occur in the synapse by neglecting opioid effect on NMDAR’s conductance (CK=0 label in Figure 4). Although, neglecting opioid’s effect on Mg2+ blocking (CMg=0 label in Figure 4) is not a critical parameter, knocking out VGCC’s activation leads to attenuation of the LTP (Icr=0 label in Figure 4). These results show that enhancement of NMDAR’s conductance besides VGCC’s activation is a fundamental parameter in inducing pathological LTP. Therefore in the next part, constants for VGCC’s activation and NMDAR’s conductance are changed for further simulations.
The left panel of Figure 5 shows the normalized value for phosphorylated CaMKII and AMPAR with a variable coefficient of NMDAR’s conductance. It can be seen that when the coefficient is 0.2, there is not any phosphorylation; however, when it is 0.3, CaMKII and AMPAR phosphorylate. It means that attenuation of LTP by manipulating NMDAR’s conductance is impossible. We can whether completely block LTP or not, and there is nothing in between. However, according to the simulation results shown in right panel of Figure 5, increment in VGCC mediated current strengthens LTP. As a matter of fact, increasing channels activation constant intensifies LTP. 
Besides opioid’s direct effect on postsynaptic neurons, change in frequency and amplitude of stimulation signal which are considered as synaptic glutamate can alter pathological LTP. The left panel of Figure 6 shows the variation of glutamate on the CaMKII and AMPAR phosphorylation process. This simulation shows that reducing glutamate amplitude by 0.01 increments from 0.2 does not lead to inhibition of LTP until glutamate amplitude reaches 0.14. At this value, phosphorylated CaMKII and AMPARs decline rapidly, and LTP does not occur. However, a gradual decrease in stimulation frequency leads to gradual descend in LTP induction shown in the right panel of Figure 6. 

4. Discussion
Our simulation results support the hypothesis that opioids consumption can lead to pathological memory formation in hippocampal neurons through induction of LTP. To explain this issue, we developed a computational model in which CaMKII and AMPAR phosphorylation processes are considered as foot prints of LTP induction. 
Experimental observations (Chen & Huang, 1992; Chen & Marine, 1991; Kow et al., 2002; Mao, 1999) indicate the opioids effects on NMDARs. Enhancement of NMDAR mediated current is a byproduct of opioids presence in synapses. Our simulation results showed that enhancement of NMDAR’s conductance is more efficient in inducing LTP compared to inhibition of Mg2+ blocking through opioids. In other words, despite the inhibitory effect of Mg2+ blocking on NMDAR mediated current, we observed that the main reason of LTP induction is enhancement of NMDA receptors conductance. Furthermore, our results suggest that prevention of pathological LTP can be achieved by blocking NMDAR’s conductance. NMDAR’s conductance acts like a none-or-all switch in inducing LTP. Therefore, it is impossible to attenuate pathological LTP by decreasing NMDAR’s conductance. However, we may completely block it using this approach. 
Bongianni et al. (1986) and Mishra et al. (2017), through their experimental observations, suggested that blocking VGCCs leads to preventing dependency on opioids. In this research, we obtained the same result but in a different way, i.e., we showed the role of VGCCs in forming pathological LTP. Our simulation results demonstrate that by inhibiting VGCCs, pathological LTP is attenuated. Also, we may infer that by decreasing VGCC activation, pathological LTP is gradually lessened. One may be able to apply this finding in cases where one intends or prefers to attenuate pathological LTP rather than blocking it. 
On the other hand, disinhibitory effect of opioid increases glutamate concentration in the synaptic cleft. Therefore, it is reasonable to assume that changes in the amplitude and the frequency of the stimulation signal in the postsynaptic neuron are due to the presence of the opioid. Simulation results demonstrate that increasing the frequency and the amplitude of the stimulation signal results in LTP induction. There are different experimental observations that support these simulation results (Akaishi et al., 2000; Capogna et al., 1993; Caudle & Chavkin, 1990; Cohen et al., 1992; Kauer & Malenka, 2007; Rezai et al., 2013). 
These observations indicate that opioid inhibits the GABAergic signal that leads to enhancement of glutamate concentration in the synapse. Furthermore, our simulation results suggest that reducing synaptic glutamate can prevent pathological LTP induction. Amplitude and frequency of glutamate release are parameters that will help in preventing or attenuating LTP. Since the source of glutamate is mainly from presynaptic neurons, inhibition of presynaptic pyramidal neurons may help prevent pathological LTP induction. On the other hand, one of the fundamental functions of the glial cells is reuptake of the synaptic glutamate. Therefore, we also suggest considering the role of glial cells in preventing pathological memory formation. This goal can be achieved through stimulating glutamate transporters. Using this approach, synaptic glutamate can be reduced, and this can result in reducing stimulation of the postsynaptic receptors. 
In this paper, we introduced a computational model for postsynaptic neurons in CA1 region of the hippocampus to examine LTP induction due to opioids. The model is consisted of a spine which produces EPSPs. NMDAR, AMPAR and VGCC’s mediated currents are also considered in the model. Calcium ions flow into the neuron through these channels and receptors. Elevation of calcium concentration leads to CaMKII and AMPAR phosphorylation. We assumed that CaMKII and AMPAR phosphorylation are foot prints of LTP induction. 
Simulation results of the presented model indicate that presence of opioids can lead to LTP induction in a situation that normal LTP does not occur. This finding supports the idea that opioids may form a memory that is directly linked to drug consumption. Therefore, one of the primary triggers of relapse in withdrawal time is this pathological memory. Furthermore, we observed that, by decreasing VGCC’s activation, NMDAR’s activation, and glutamate concentration in the synapse, pathological LTP reduces and it is even possible to prevent the pathological LTP. 

Ethical Considerations
Compliance with ethical guidelines
Since the research was a simulation one, the authors did not need to get any ethical approvals.

Funding
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Conflict of interest
The authors declared no conflict of interest.


References