Mathematical Modelling MODELING OF AMBIENT GLUTAMATE CONCENTRATION MEASUREMENT IN THE MAMMALIAN NERVOUS SYSTEM D. Shchepakin1 , M. Kavanaugh2 , L. Kalachev1 1 Department of Mathematical Sciences, University of Montana, Missoula, MT, USA 2 Center for Structural and Functional Neuroscience, University of Montana, Missoula, MT, USA Abstract. A neuron is an electrically excitable cell that processes and transmits information through electrical and chemical signals. Neurons connect and pass signals to other cells through the struc- ture called synapse. We focus on synapses through which the signals are transferred by signaling molecules called neurotransmitters. One of the predominant excitatory neurotransmitters in the central ner- vous system of the mammals, including humans, is glutamate. It is directly or indirectly involved in most brain functions. However, the excessive stimulation of the glutamate receptors is toxic to neurons, therefore it is important to rapidly clear the glutamate from the extra- cellular space and keep its concentration low. Glutamate transporters play a crucial role in regulating glutamate concentration in synaptic clefts. Thus, it is important to understand the mechanisms underlying this process. We describe measurement of the glutamate concentra- tion in the extracellular space. It is important to estimate the baseline glutamate concentration to use it in future models and studies. How- ever, two existing methods of measuring the glutamate concentration in the extracellular space give inconsistent results with about 100 fold difference. We construct the model of the process of the glutamate concentration measurement in order to explain that discrepancy. Keywords: glutamate transport, ambient neurotransmitter, tonic signaling, microdialysis Citation: Shchepakin D, Kavanaugh M, Kalachev L. Modeling of am- bient glutamate concentration measurement in the mammalian ner- vous system. CEUR Workshop Proceedings, 2016; 1638: 717-730. DOI: 10.18287/1613-0073-2016-1638-717-730 Introduction A neuron is an electrically excitable cell that processes and transmits informa- tion through electrical and chemical signals. A typical neuron has a cell body, multiple dendrites and one axon (Figure 1). Because neurons are electrically excitable, change in a cross-membrane potential could cause the activity of volt- age dependent ion channels. If this change is high enough, it will cause a one way electrochemical pulse along the axon, called action potential. When action Information Technology and Nanotechnology (ITNT-2016) 717 Mathematical Modelling Shchepakin D, Kavanaugh M, Kalachev L. Modeling... Fig. 1. Neuron and synapse structures (http://commons.wikimedia.org/wiki/File:Chemical synapse schema cropped.jpg) potential reaches the end of the axon it will activate the series of reactions, which transfer the signal to the adjacent cell the axon is connected to. This connection is called the synapse connection and the connection structure is called a synapse. The cell that initiates the signal is called a presynaptic cell, and the one that receives the signal is called a postsynaptic cell. Depending on the type of reac- tions that process in the synapse to transfer the signal, there are two different types of synapses: chemical synapses and electrical synapses. We will focus on the chemical synapses. The end of the axon in the synapse connection contains vesicles with signaling molecules called neurotransmitters, which are released into the thin space, called synaptic cleft, between presynaptic and postsynaptic cells when action potential hits. The postsynaptic cell possesses the receptors that are triggered by the neurotransmitters, which allows the postsynaptic cell to that start the cascade of the reactions leading to the further activity. We will consider the dynamic of the glutamate, which is the predominant excita- tory neurotransmitter in the central nervous system of the mammals, including humans. It is directly or indirectly involved in most brain functions. However, Information Technology and Nanotechnology (ITNT-2016) 718 Mathematical Modelling Shchepakin D, Kavanaugh M, Kalachev L. Modeling... the excessive stimulation of the glutamate receptors is toxic to neurons, there- fore it is important to rapidly clear the glutamate from the extracellular space and keep it low. Glutamate transporters play a crucial role in regulating gluta- mate concentration in synaptic clefts. Thus, it is important to understand the mechanism of this process. Methods of measuring the glutamate concentration in a brain The knowledge of the concentration of glutamate in the extracellular space is cru- cial for applying theoretical models describing its dynamic. Two different meth- ods exist to estimate this value: microdialysis and electrophysiological methods. However, microdialysis measurement turns out to be about 100-fold higher than the electrophysiological measurement [1]–[4]. The reasons for this discrepancy are yet unknown. Fig. 2. Schematic illustration of a microdialysis probe (http://commons.wikimedia.org/wiki/File: Schematic illustration of a microdialysis probe.png) Microdialysis is an invasive sampling technique that is widely used to measure the concentrations of free analyte in the extracellular fluid. The idea is to insert a microdialysis catheter with semipermeable walls into the tissue (Figure 2). Be- cause of the perfusion, one is able to measure the concentration of the substance of interest after some time when the stable state of the system has been reached. It was suggested that because of the invasiveness of the microdialysis method the glutamate transporters activity could be reduced in a small region adjacent to the dialysis probe, which might lead to the overestimation of the glutamate concentration in the extracellular space. We will construct a mathematical model of this process under several assump- tions in order to investigate the dynamic and the steady state of the glutamate concentration inside the dialysis probe and in a small region nearby. Our goal is to find a possible explanation for the discrepancy of the two mentioned above methods using the constructed model. Information Technology and Nanotechnology (ITNT-2016) 719 Mathematical Modelling Shchepakin D, Kavanaugh M, Kalachev L. Modeling... Mathematical model Model assumptions Before constructing our model we need to formulate all assumptions on which it will be based. It is reasonable to assume that glutamate molecules do not appear from or dis- appear to nowhere. Thus, all changes in glutamate concentration must follow the conservation law. Moreover, because of the nature of the process, one can expect the model may be formulated in the form of a diffusion partial differen- tial equation. Therefore, the first Fick’s law can apply. Furthermore, for the simplicity of the model we will assume that the diffusion coefficient in this law is constant. Now, we discuss some basic assumptions on cells and on glutamate transporters. To transport any glutamate into the cell, a glutamate molecule must bind to a transporter. After that the glutamate molecule can either unbind from the transporter returning the system to the previous state or penetrate inside the cell via transporter, again unbinding from the transporter. There is an evidence that binding and unbinding rates of glutamate molecules and glutamate transporters are much faster than the transfer rate (also called turnover rate) of glutamate molecules into the cell by the transporters [5]. Additionally, a leak of glutamate from the cells into the extracellular space by simple diffusion is reported [6]. The evidence exists that the concentration of the glutamate is much higher inside the cells than in the extracellular space [7], which suggests that the leak could be assumed to be constant. Next, let us put several assumptions about the microdialysis catheter and the region of the treatment. First, we will assume that the dialysis probe is a perfect circular cylinder. Moreover, there are no sources or sinks of the glutamate inside the probe. We assume that only the walls of the catheter are permeable for the glutamate, that is, there are zero fluxes of glutamate through the top and the bottom of the cylinder. We will also assume that the conditions inside and outside the probe are symmetric with respect to the tube axis and along it. Finally, we will assume that the invasiveness of the treatment caused the re- duction in transporting rate of the transporters. The closer the cells are to the location of the insertion, the bigger the reduction in transporting rate. We will assume, that the cumulative reduction of the transporting rate of all transporters can be approximated by a Gaussian function as a function of a distance from the location of treatment. Model equation derivation One can imagine that modeling of the dynamics of the glutamate concentra- tion can be represented as the movement of many individual molecules in some volume. Let us call this region Ω, the open subset of Rn , where n ≥ 1. Of course, the most interest for us will be the case n = 3. Moreover, because of the assumed symmetry of the Ω it is possible to map it on the subset of a space of lower dimension, but, for now, we will consider the case of n = 3. Let us introduce more notations, that we will use to construct our mathematical model. • Let V be an arbitrary compact volume such that V ⊂ Ω and ∂V be a piecewise smooth boundary of the region V . Information Technology and Nanotechnology (ITNT-2016) 720 Mathematical Modelling Shchepakin D, Kavanaugh M, Kalachev L. Modeling... • The density function u(t, x), where x ∈ Ω is a point in 3-dimensional space and t ≥ 0 is time. The units of the values of this function correspond to the volume density, that is, the amount of the glutamate per unit volume. We will assume, that the density function u(t, x) is sufficiently smooth. • Let φ(t, x) be the flux of the glutamate at the location x ∈ Ω, at time t ≥ 0. It measures the amount of the glutamate crossing the point x, at time t ≥ 0. The units of the values of this function are the amount of glutamate per unit area, per unit time. One should notice, that φ : R × R3 → R3 , that is, the flux function is a vector function. • Let f (t, x) be the rate at which the glutamate appears or disappears per unit volume at the location x ∈ Ω, at time t ≥ 0. We will call the function f a source if it positive, sink if it is negative and source-sink if it takes both positive and negative values or if the sign of its values is unclear. Let’s now consider the dynamics of a glutamate amount within the arbitrary volume V . According to the conservation law, the rate of change of the amount of the glutamate in that volume must be equal to the rate at which the glutamate amount flows in through the boundary ∂V minus the rate at which it flows out through the boundary ∂V plus the rate at which it appears within the volume V minus the rate at which it disappears within the volume V . Therefore, one can write this in the mathematical form, using the introduced notations, as follows Z Z Z d u(t, x)dv = − φ(t, x)da + f (t, x)dv. (1) dt V ∂V V Now, applying the Divergence theorem to the first term on the right hand side of the equation (1), we obtain Z Z Z d u(t, x)dv = − ∇ · φ(t, x)dv + f (t, x)dv. dt V V V Using Leibniz integral rule, we get Z Z Z ∂ u(t, x)dv = − ∇ · φ(t, x)dv + f (t, x)dv. ∂t V V V Next, moving all term to the left hand side and combining all the integrals, we can write Z   ∂ u(t, x) + ∇ · φ(t, x) − f (t, x) dv = 0. (2) ∂t V Because the equation (2) must hold for any volume V , and because the functions under the integral are continuous, we arrive at ∂ u(t, x) = −∇ · φ(t, x) + f (t, x). (3) ∂t Fick’s first law, which states that the flux goes from regions of high concentration to regions of low concentration with a magnitude that is proportional to the concentration gradient can be written in mathematical form as φ(t, x) = −D(t, x)∇u(t, x), Information Technology and Nanotechnology (ITNT-2016) 721 Mathematical Modelling Shchepakin D, Kavanaugh M, Kalachev L. Modeling... where D is a diffusion coefficient. We will assume, that the diffusion coefficient doesn’t depend on time and location, that is, D is a constant. Thus, substituting the flux expression into the equation (3), we obtain ∂ u(t, x) = D4u(t, x) + f (t, x). (4) ∂t Our next step is to derive the expression for the function f . We can represent this function as a combination of two processes: one process is related to the chemical kinetics of the glutamate molecules because of the reacting with gluta- mate transporters and another is related to the leak of the glutamate molecules from the cells into the extracellular space. We write f (t, x) = fck (t, x) + fleak (t, x). (5) Because of the assumptions we made, that leak is invariable in space and time: fleak ≡ KL , where KL is a constant. Therefore, the equation (5) can be written as f (t, x) = fck (t, x) + KL . (6) Let us now consider time dependent chemical kinetics of glutamate molecules and transporters in an arbitrary volume W ⊂ Ω. To transfer a glutamate molecule inside the cell, the molecule must bind to the glutamate transporter first to form a complex. Moreover, we assumed that the bound glutamate molecule could unbind from and free the transporter before it will transfer the molecule into a cell. The corresponding chemical kinetics scheme for this process can be written as follows k+ kV Gout + T r − − )* −− C, C −−→ T r + Gin , k− where Gout is extracellular glutamate, Gin is intercellular glutamate, T r is a transporter and C is a glutamate-transporter complex, k+ and k− are the binding and unbinding rate constants respectively and kV is the rate of transfer of the glutamate into a cell (thus, the glutamate goes into the cell and complex becomes a free transporter). Applying the Law of Mass Action, we can write the system of the differential equations that describes those reactions as follows: d[Gout ] = −k+ [Gout ][T r] + k− [C], dt d[T r] = −k+ [Gout ][T r] + (k− + V )[C], dt (7) d[C] = k+ [Gout ][T r] − (k− + V )[C], dt d[Gin ] = V [C], dt where [A] is a notation for the concentration of substance A and d[A] dt is the time derivative of the concentration of substance A. Information Technology and Nanotechnology (ITNT-2016) 722 Mathematical Modelling Shchepakin D, Kavanaugh M, Kalachev L. Modeling... We will denote the initial conditions as [Gout ](0) = G∗out , [T r](0) = T r∗ , [C](0) = C ∗ , [Gin ](0) = G∗in . One should notice that the total number of the glutamate transporters and, therefore, the concentration of the transporters, which is the sum of concen- tration of free transporters and bound transporters (complexes), should remain constant. Indeed, d[T r] d[C] + = 0, dt dt so [T r] + [C] = const = [T r](0) + [C](0) = T r∗ + C ∗ = N, then [T r] = N − [C]. (8) Now, according to our assumption, the binding/unbinding rates of the glutamate molecules to transporters are much faster than the rate at which the glutamate molecules are transferred inside the cells (via the glutamate transporters). We can write this assumption in the mathematical form using our notations to get k+ N ∼ k−  kV . Let’s introduce a small parameter 0 < ε  1 and a new notations kf + = εk+ , k− = εk− . f Then we have + N ∼ k− ∼ kV . kf f Using these new notations, we can rewrite the system (7) in the following form d[Gout ] ε = −kf + [Gout ][T r] + kf − [C], dt d[T r] ε = −kf + [Gout ][T r] + (kf − + εV )[C], dt (9) d[C] ε = kf + [Gout ][T r] − (kf − + εV )[C], dt d[Gin ] = kV [C]. dt Our next step is to apply Thikhonov’s and Vasil’eva’s theorems [8] to this per- turbed problem to obtain the reduced version of the original system (7). Because we are not really interested in the behavior of the system near initial instant of time, we will consider only the, so called, regular part of the asymptotic ex- pansion. For the sake of simplicity, we will consider only leading order approx- imation. Thus, we can represent the solution of the system (9) and, therefore, Information Technology and Nanotechnology (ITNT-2016) 723 Mathematical Modelling Shchepakin D, Kavanaugh M, Kalachev L. Modeling... original system (7) as [Gout ] = [Gout ] + O(ε), [T r] = [T r] + O(ε), [C] = [C] + O(ε), [Gin ] = [Gin ] + O(ε). Substituting these expressions into the system (9) and considering only equations of zero order approximation ε, we will have 0 = −kf + [Gout ] [T r] + kf − [C], 0 = −kf + [Gout ] [T r] + kf − [C], 0 = kf + [Gout ] [T r] − kf − [C], ˙ [Gin ] = kV [C], which leads to kf + [Gout ] [T r] = kf − [C], ˙ [Gin ] = kV [C]. After substituting (8) into the first equation above, we write + [Gout ] (N − [C]) = k− [C], kf f and + N [Gout ] kf [C] = . − + k+ [Gout ] kf f Thus, the reduced system is N [Gout ] [C] = , Kd + [Gout ] (10) ˙ [Gout ] [Gin ] = N kV , Kd + [Gout ] k− where Kd = . k+ The second equation of the system (10) describes the leading order approxima- tion of the glutamate flow rate from the extracellular space into the cells due to the activity of the glutamate transporters. It depends on the concentration of transporters N , their turnover rate kV , the binding-unbinding rate constants k± and the concentration of the glutamate outside the cells. Because of the Information Technology and Nanotechnology (ITNT-2016) 724 Mathematical Modelling Shchepakin D, Kavanaugh M, Kalachev L. Modeling... assumptions we made about the resulting transfer rate of all transporters we can now express function fck (t, x) in the notations of the equation (4) as u(t, x) fck (t, x) = J(l) , Kd + u(t, x) where l is a distance between the location x and the axis of the dialysis probe. Also,    0, 0 ≤ l ≤ L,       J(l) = (d − L)2 (11) − 2σ 2 J 1 − e , l > L,    max        where Jmax is a cumulative transfer rate coefficient in a healthy tissue, L is a radius of the microdialysis catheter, and σ is some parameter that corresponds to the amount of influence of the treatment on the turnover rate of the transporters. Thus, the equation (4) can be written as ∂ u(t, x) u(t, x) = D4u(t, x) − J(d) + KL . (12) ∂t Kd + u(t, x) Because of the assumed symmetry and non-permeability of the top and the bot- tom of the cylindrical dialysis probe, the equation (12) can be greatly simplified and reduced to a lower dimension equation. Fig. 3. The cylindrical probe. It is more convenient to consider the equation (12) in cylindrical coordinates. So, we need to change variables using the following expressions, assuming that x = (x, y, z) x = r cos φ, y = r sin φ, Information Technology and Nanotechnology (ITNT-2016) 725 Mathematical Modelling Shchepakin D, Kavanaugh M, Kalachev L. Modeling... where r is the Euclidean distance from the z axis to the point x and φ is the azimuth of the point, that is, the angle between the reference direction on the xy plane and the line from the origin to the projection of point x on the plane. Moreover, we make the origin of the coordinate system be in the center of a lower base of the cylindrical probe with z axis being coinciding to the axis of the cylinder (see Figure 3). The equation (12) will transform into ∂ u(t, r, φ, z) u(t, r, φ, z) = D4u(t, r, φ, z) − J(r) + KL , (13) ∂t Kd + u(t, r, φ, z) where 4u(t, r, φ, z) is a Laplace operator in a cylindrical coordinates acting on u: 1 ∂2u ∂2u   1 ∂ ∂u 4u(t, r, φ, z) = r + 2 2 + 2. r ∂r ∂r r ∂φ ∂z Because of the assumptions we made about symmetry with respect to the axis of the cylindrical probe and along that axis, the function describing glutamate concentration will not depend on φ and z, that is, u(t, r, φ, z) = u(t, r). Therefore, the Laplace operator in our case will be just ∂ 2 u 1 ∂u   1 ∂ ∂u 4u(t, r) = r = + . r ∂r ∂r ∂r r ∂r Now we can rewrite equation (13) as ∂ 2 u 1 ∂u   ∂ u u=D + − J(r) + KL , (14) ∂t ∂r r ∂r Kd + u where u = u(t, r). Finally, we need to derive boundary and initial conditions for the constructed model equation. Again, using the symmetry assumptions, we can state that there is no flux through the center of the cylindrical catheter. We can write this as ∂u = 0. (15) ∂r r=0 Next, let us denote the concentration of glutamate in the extracellular space by us . We expect the microdialysis treatment to affect the concentration of glutamate only in some small region, therefore, u(t, ∞) = us . (16) The initial conditions for the concentration of the glutamate: some specified concentration u∗ inside the dialysis tube, which can be zero or not, and us outside the tube:  ∗  u , 0 ≤ r ≤ L, u(0, r) = (17) us , r > L.  Information Technology and Nanotechnology (ITNT-2016) 726 Mathematical Modelling Shchepakin D, Kavanaugh M, Kalachev L. Modeling... Let us write down the final version of constructed model (14), (11), (15), (16), (17):  2  ∂u ∂ u 1 ∂u u =D + − J(r) + KL , ∂t ∂r r ∂r Kd + u    0, 0 ≤ r ≤ L,       J(r) = (r − L)2 (18) −     Jmax 1 − e  2σ 2 , r > L,   ∂u = 0, u(t, ∞) = us , ∂r r=0  ∗  u , 0 ≤ r ≤ L, u(0, r) = us , r > L.  This problem cannot be solved analytically because of the nonlinear term in the right hand side of the diffusion equation. Next, we will consider the numerical method that will allow us to solve this problem computationally. Implementation and results Let us introduce physically realistic constants that will be substituted into the model equation (18) [9]. L = 1500µm, σ = 200µm, us = 25nM, Kd = 25µM, Jmax = 2.5mM sec−1 , KL = 2.5µM sec−1 , D = 700µm2 sec−1 . We will use different initial concentrations of glutamate inside the cylindrical probe (u∗ ) to find the cylindrical volume around the probe, so that the radius of the base of the cylindrical volume R is as small as possible, but large enough, so that it ”can still play” the role of an ”infinity” in our model. We can say that the value R1 of R is large enough if substituting some R2 > R1 into the model will not change the values of the solution u(t, r), provided by the numerical method, on the interval [0, R1 ]. We should notice that we assumed the ”true” concentration of glutamate in the extracellular space to be 25nM . This values is suggested by the electrophysiology method of measurement of glutamate concentration. We ran the simulation for different choices of values of R and u∗ and obtained the results, shown in Figures 4, 5, and 6. Information Technology and Nanotechnology (ITNT-2016) 727 Mathematical Modelling Shchepakin D, Kavanaugh M, Kalachev L. Modeling... Fig. 4. Spatial profiles of glutamate concentration at the initial time t = 0 seconds, at t = 5400 seconds and at the steady state. The initial concentration in the tube is u∗ = 10mM . First row corresponds to the choice of R = 3000µm, second row corresponds to the choice of R = 6000µm. Fig. 5. Spatial profiles of glutamate concentration at the initial time t = 0 seconds, at t = 5400 seconds and at the steady state. The initial concentration in the tube is u∗ = 0M . First row corresponds to the choice of R = 3000µm, second row corresponds to the choice of R = 6000µm. Information Technology and Nanotechnology (ITNT-2016) 728 Mathematical Modelling Shchepakin D, Kavanaugh M, Kalachev L. Modeling... Fig. 6. Function of glutamate concentration at the center of the probe as a function of time. The initial concentration in the tube was u∗ = 10mM. We can notice in Figures 4 and 5 that the resulting concentration looks be the same for R = 3000µm and R = 6000µm. Computations shown no significant difference between the solutions for those two values of R. Therefore, R = 3000µm can be used to represent ”infinity” in our numerical model. Moreover, we note that no matter what initial conditions are, the system always reaches the same steady state. The only difference is the time interval within which this steady state is reached. Therefore, one could use any initial conditions to analyze the steady state of the system, given enough time passed. As we can see on the graphs, the concentration of glutamate at the steady state inside the microdialysis probe is constant, i.e. the glutamate there is uniformly distributed. In Figure 6 the time dependence of glutamate concentration at the center of the cylindrical catheter (r = 0) is shown. Let us now compare the values of concentration inside and outside of the dialy- sis catheter at the steady state. The concentration of glutamate outside of the catheter at ”infinity” is just the glutamate concentration in the extracellular space us . By our assumption, it’s the concentration provided by the electro- physiology method, and it equals to 25nM . Under that assumption the model predicts about 3.5µM for the glutamate concentration inside the catheter, which is within the range of the usual values estimated via microdialysis method. To sum up, we constructed a model of the microdialysis method under sev- eral assumptions. Assuming that the true concentration of the glutamate in the extracellular space is provided by the electrophysiology method, the model predicts the values of the microdialysis measurements that match the values ob- tained by the real measurements. Therefore, we provided possible explanation for the discrepancy in measurements between microbialysis and electrophysiol- ogy methods. Information Technology and Nanotechnology (ITNT-2016) 729 Mathematical Modelling Shchepakin D, Kavanaugh M, Kalachev L. Modeling... References 1. Zerangue N, Kavanaugh MP. Flux coupling in a neuronal glutamate transporter. Nature, 1996; 383: 634-637. 2. Levy LM, Warr O, Attwell D. Stoichiometry of the glial glutamate transporter GLT-1 expressed inducibly in a Chinese hamster ovary cell line selected for low endogenous Na+-dependent glutamate uptake. J. Neurosci, 1998; 18: 9620-9628. 3. Benveniste H, Drejer J, Schousboe A, Diemer NH. Elevation of the extracellular concen- trations of glutamate and aspartate in rat hippocampus during transient cerebral ischemia monitored by intracerebral microdialysis. J. Neurochem, 1984; 43(5): 1369-1374. 4. Lerma J, Herranz AS, Herreras O, Abraira V, Martn del Ro R. In vivo determination of extracellular concentration of amino acids in the rat hippocampus. A method based on brain dialysis and computerized analysis. Brain Res., 1986; 384(1): 145-155. 5. Wadiche JI, Arriza JL, Amara SG, Kavanaugh MP. Kinetics of a human glutamate trans- porter. Neuron, 1995; 14: 1019-1027. 6. Lehre KP, Danbolt NC. The number of glutamate transporter subtype molecules at glutamatergic synapses: chemical and stereological quantification in young adult rat brain. J. Neurosci., 1998; 18(21): 8751-8757. 7. Niels C. Danbolt. Progress in Neurobiology, 2001; 65(1): 1-105. 8. Tikhonov AN, Vasil’eva AB, Sveshnikov AG. Differential Equations. Springer Series in Soviet Mathematics. Berlin, Heidelberg, New York, Tokyo: Springer-Verlag, 1985. 9. Wadiche JI, Arriza JL, Amara SG, Kavanaugh MP. Kinetics of a human glutamate trans- porter. Neuron,1995; 14(5):1019-1027. Information Technology and Nanotechnology (ITNT-2016) 730