CC BY
45
ЭЛЕКТРОНИКА
UDC 517.958:536.71
MATHEMATICAL MODELS AND NUMERICAL SIMULATION OF DIFFUSION-REACTION PROBLEMS OF BRAIN CHEMISTRY
OLEINICK A.I., AMATORE C.A., SVIR I.B. * 1
The processes of release of NO by a neuron, its diffusion through the brain media and collection by the electrode are numerically simulated using the conformal mapping approach initially proposed in the article (C. Amatore, A.I. Oleinick, I.B. Svir//J. Electroanalytical Chemistry. V. 553, 2003, pp. 49-61). Fitting simulated currents to experimental ones yielded the parameters of the neuron activation function and the value of the surface concentration of NO.
1. Introduction
The intimate integrated mechanisms which control the local brain and cerebellum activities are challenges in modern biochemistry and they have been within the scope of biochemical community for the last several decades. Understanding of such basic mechanisms is very important and provides new insights in how the brain works and therefore will have a lot of applications, especially in medicine, biochemistry etc.
In this connection, Functional Imaging Magnetic Resonance (IRMf) provides important clues on how the brain performs its main function. It is the great advantage of IRMf to monitor efficiently and in-time in a perfectly non-invasive way the differential status (ca. one percent at most) of blood circulation flows in the living brain. This allows 3D-identification of the precise zones which experience larger blood flows without any significant interference.
For many medical, therapeutic, biological and behavioural investigations the IRMf data are used within the framework of serendipitous empirical correlations involving many other biological and instrumental factors. Therefore, it suffices for most of the present IRMf studies that the existence of overall correlations is established between the studied phenomena and any variations in the IRMf 3 D-mapping of the living brain, though the fundamental biological and physicochemical reasons giving rise to these correlations remain fundamentally occulted, since there is still no understanding of the mechanism(s) on which IRMf analyses are crucially rooted. The mechanisms which provoke a local increased blood supply, viz., a dilatation of blood capillary vessels, in the active areas or their kinetics remain unknown [1-3].
In this paper, the process of nitric oxide release from an active neuron is studied, which is believed to control the dilatation of blood capillary vessels and hence the
blood supply to active areas of the living brain. The NO released by a single patch-clamped neuron is monitored using an ultramicroelectrode introduced into a slice of living rat’s brain.
We used experimental results of the detection of nitric oxide during neuronal activity obtained recently via electrochemical methods [4]. In order to extract information from these experimental data and in order to confirm the existence of the interrelationship between the release of nitric oxide by the active neuron and dilatation of the blood vessels the underpinning theory should be developed. The aim of this paper is to elaborate on such theory, viz. on the corresponding mathematical models and an efficient approach to their numerical solution.
2. Theory
2.1. Mathematical models
Consider a brain slice with a neuron within it and a blood vessel passing through the slice as shown in Fig. 1a, i.e. the generatrices of a cylinder representing the blood vessel are parallel to the y-axis. Fig. 2a shows a section of the brain slice taken in the middle of its height, i.e. it coincides with the (x, y) plane in Fig. 1a. The neuron is of radius r and the line of intersection of its surface with the secant plane (x,y) is depicted as Ц in Fig. 2a. The line of intersection of the cylindrical surface of the blood vessel with the secant plane is a straight line depicted as Г2 in Fig. 2a. The neuron and blood vessel are separated by a gap of the size g .
Fig. 1. Schematic representation of the brain slice: a — initial view of the system; b — layout of the system after insertion of the ultramicroelectrode into the slice
As has been mentioned above active neurons release some amount of NO, which reaches its target (i.e. the vessel) through diffusion within the brain slice media (designated with grey colour in Fig. 1). When molecules of NO reach the blood vessel they are absorbed at the vessel surface and cause the dilatation of the latter.
18
РИ, 2005, № 3
However, when an electrode is introduced into the brain slice to monitor the neuron’s activity (see Fig. 1b) most of NO is intercepted before it can arrive to the vessel surface. Therefore in our model we do not take into account the presence of the vessel and limit ourselves with exclusive consideration of the neuron and the electrode (as depicted in the cross-section in Fig. 2b).
It is also established that released NO interacts with the oxygen present in the brain slice media and this process may be described by the following rate law [5, 6]:
dc
dt
k2c2 ,
(1)
where c0 is the maximum concentration at the neuron surface; n is the unit normal vector to the electrode surface; k(t) is the neuron activity function, Sym — represents the part of the boundary corresponding to the symmetry axis,
Sym = {x:xє (-да, -q)u(rbq + g)u(p + Г2, + да)},
p = q + g + Г2 . The Dirichlet boundary condition is imposed at the surface of the electrode since the electrode potential was chosen to be at the plateau of the NO oxidation wave.
where c is the concentration of NO; k2 is the second order rate constant, k2 = 9.3 x 106 molcm_3 s_1.
The thickness of the brain slice in the z-direction is much smaller than its dimensions in x - and y -directions, so we assume that the diffusion of NO is essentially two-dimensional and neglect the mass transport of NO in the z -direction. This means that simulations can be performed for the cross-sections shown in Fig. 2a,b.
The above-mentioned assumption allows considering such mathematical model as a first level approximation. However, we believe that such simplification of the mathematical model is useful since it allows building the framework for the problem at hand and revealing principal properties and dependences.
Thus, the mathematical model of the mass transport of NO in the neuron-electrode system can be written as following:
= D d t
d 2c
d 2c
dx2 d y2
- k2 c
2
(2)
where D is the NO diffusion coefficient and t is the time. The corresponding initial and boundary conditions are:
t = 0 : V x,y; c = 0 ;
^ c
t > 0 : Vx, yє Г1; D— = -k(t)[c0 -c] (neuron) d n
3 c
y = 0 , x є Sym ; — = 0 (symmetry axis)
d y
Vx, y є Г2; c = 0; (electrode)
x2 + y2 ^ да; c = 0 , (infinity)
(3)
Considering of the pure neuron-vessel system allows one to calculate the amount of NO absorbed by the vessel (on the basis of the parameters obtained from fitting experimental data for the neuron-electrode system) and reveal the dependence between this amount and the dilatation of the blood vessel. However this case will not be considered here and is the subject for the subsequent contribution.
Herein we wish to emphasize that both of the mentioned systems can be treated in a similar way due to the application of the conformal mapping (see below). The mathematical model for the neuron-vessel case is essentially the same except for the symmetry axis which
becomes Symves = { x :x є (-да, - q) u (q, q + g)} and the boundary condition for electrode in (3) which is replaced by the following boundary condition for the vessel:
3 c
Vx, ye Г2 ; D— = kabs c , (vessel) (4) О n
where kabs is the rate constant ofabsorption of NO by the surface of the blood vessel. The cross-section of the neuron-vessel system is depicted in Fig. 2a.
2.2. Conformal mapping
In order to generalise the analysis as well as the results of our simulations we introduce the following dimensionless variables and parameters:
_ c x y D t
c0 r1 r1 r12
KM = ^; K2 = k2 c0r12 .
Kabs _
D
kabs r1
D
D
R1 = -1 = 1; G ;
£L
r1
_g
r1
R2 = 2Г , P = -P = 1 + G + R2
(5)
Fig. 2. Cross-sections of the brain slice along the (x,y) -plane which correspond to the situations depicted in Fig.1a,b respectively
РИ, 2005, № 3
19
The domain shown in Fig. 2b can be mapped onto a closed rectangle (see Fig. 3) by the following conformal map [7]:
ae (a-b)cosp + be^
X —-----------------------
2(ch^ - cos p)
Y = (a - b)sin p 2(ch^ - cos p) ’
where
(6)
(1+ P2 - r2) -^(1+ P2 - R2)2 - 4P 2 b i;
a =--------------------------------, b = 1/a
2P
are the parameters defined by the geometry of the system.
The non-uniform grid shown in Fig. 2b corresponds by the conformal map (6) to the uniform rectangular grid shown in Fig. 3.
The transform (6) also maps the domain depicted in Fig .2a onto a rectangular area since the vessel boundary,
Г2, can be considered as a circle of infinite radius with the centre located at infinity. Thus, the vessel boundary is mapped onto the segment (£,2 = 0 , л<р< 2n ) in the conformal space, while the left boundary (§1) remains unchanged. However, it should be noted that for obtaining correct values of the parameters a, b of the mapping (6), for the neuron-vessel case, they should be calculated as if there existed another circle of the same radius, Ri, as a neuron with the centre located at X = 2 + G.
Fig. 3. Conformal space with rectangular grid corresponding to curvilinear grid in Fig. 2b
The non-uniform grids shown in Fig. 2 correspond by the conformal map (6) to the uniform rectangular grid shown in Fig. 3.
2.3. Mathematical model in conformal space
The conformal mapping (6) was applied to the mathematical model (2-3) after its normalisation using the dimensionless variables defined in (5). The model has the following view in the conformal space:
dC_D* d^C | d2C дт |_d§2 dp2
where D* = e-2^ (l + e2^ - 2e^ cosp)2/(a - b)2 is the Jacobian of the transformation (6).
20
- K2C
2
(7)
The initial and boundary conditions for the neuron-electrode case are:
T = 0: V5, P ; C = 0 ;
x> 0 : | = ^ , л<р< 2n ,
d C
h§r| = _K(x) t1 _ C] (neuron)
d C
§1 < § < § 2 , h = n, — = 0 (symmetry axis)
dp
d C
|1 2 (0), p = 2n, — = 0 (symmetry axis)
dp
'% = '%2 , ^<p<2n , C = 0 (electrode)
| = 0 , p = 2% , C = 0 , (infinity) (8)
where h^ = 2 (cosh\ - cos p)/(b - a).
The simulation area for the neuron-vessel case (as mentioned above) changes in the £ -direction to §1 <^< 0 and the boundary condition (4) rewrites as:
3 c
£ = 0, rc < P < 2% ; hл — = Kabs C . (vessel) (9)
oq
2.4. Neuron activity function
Experiments show that neuron response resembles a ‘hat’ function [4]. This prompts the idea that the neuron activity function has a very similar shape. We use the following function to describe the changes of rate constant which control the flux at the neuron body surface (see Fig. 4):
k(t) = k0 exp[a(t~^ x exp[-P(t-12)! (10)
01 + exp [a (t -11)] 1 + exp[-p (t -12)], ( )
where k0 set the required amplitude for the rate
constant; t1 and t2 are the times corresponding to the “halfwaves” of the hat function (indicated with dashed lines in Fig. 4); the coefficients a and p determine the rates of function change in the vicinity of the points 11
and t2 respectively. The dimensionless form of neuron activity function used in (8) is introduced in (5).
Fig. 4. A schematic view of the neuron activity function
РИ, 2005, № 3
2.5. Current calculations
The general expression for the evaluation of the current at an electrode is:
i = nFD f-^dS
S5n =
(11)
where S is the electrode surface, f is the Faraday constant, n is the number of transferred electrons.
For our model (Fig. 2b) the equation (11) will rewrite as follows:
3 c
i = 2nFDhel j —^2
Г2 3 n
(12)
where hei is the height of the electrode within the brain
slice as depicted in Fig. 5 and Г2 is the boundary of the electrode in considered cross-section (Fig. 2b), factor 2 in eq. (12) is due to the symmetry of the system.
The eq. (12) in conformal space rewrites in the following view:
2к d C
i = 2nFD cohei J —dp. л
(13)
Fig. 5. Geometry of the neuron-ultramicroelectrode system in a cross-section transversal to the brain slice
3. Simulations
Numerical solution of the problems (7)-(8) was performed with the ADI method, which showed a great efficiency coupling with conformal mapping technique [8-10] for the solution of 2D electrochemical problems at various symmetrical microelectrode geometries. Simulation results quoted in this paper were generated by programs written by the authors in Borland Delphi 7 Enterprise Edition and executed on a PC equipped with Intel Pentium IV processor and 512 Mb of RAM.
Despite the fact that the duration of an experiment is extremely long (ca. several hundred seconds) the uniform time grid can be used. This can be explained by the absence of initial singularity (which consists in abrupt changes of flux or concentration at the active surface at t = 0 and is the reason of instabilities in numerical solution) and the fact that concentration changes occur rather slowly — a significant change takes more than several tenths of a second.
In some cases experimental signals observed at the electrode reveal non-negligible fluxes at times that are sufficiently close to the initial moment (t = 0), which means that the flux at the neuron is also non-negligible. The latter causes the initial time singularity. In this case experimental data may be shifted towards more positive
РИ, 2005, № 3
values oft as well as the neuron activity function until the latter provides a negligible flux at the neuron surface at initial times. This can be achieved due to virtually exponential rate of decay of the activity function given
in eq.(10) for t < tj and t > t2 .
4. Results and discussions
The experimental datasets were fitted with the following physicochemical and computational parameters: D = 10_5 cm2 s_1 , hel = 5pm , N^x Nqx Nx =
120 x 120 x 10 000 . The maximum concentration at the neuron surface (c0), the amplitude ofthe time depended rate constant (k0) and the shape ofthe activity function (which is controlled by the parameters a, P , F , t2) were varied during the fitting. The geometry of the system was set as required in each case.
The numerical results given in the next three figures (Figs.6-8) are presented in the following way. The “noisy” curves are experimental data; smooth solid curves are obtained by the numerical solution of the model described by eqs. (7)-(8), (13) with the parameter values indicated precisely in each case. The curves designated with symbols correspond to numerically
simulated currents (eq. (13)) obtained with K2 = 0 in eq. (7), i.e. without taking into account the influence of homogeneous interaction between NO and oxygen. The latter curves were obtained with the same set of best fit parameter values which were obtained in the case when the homogeneous kinetics was included. This was done in order to clarify the contribution of the homogeneous kinetics.
Fig. 6. Fitting results for the system with the dimensions: q = 5 pm , g = 10 pm and Г2 = 10 pm
Fig. 6 demonstrates the fitting results for the system with the following geometrical dimensions: q = 5 pm , g = 10 pm and r2 = 10 pm . The best fit parameters are the following: c0 = 6.2x10_1° molcm-3 ,
k0 = 0.1cms_1, t1 = 146 s , t2 = 208 s , a = 0.08 , P = 0.4 . Similarly, Fig. 7 shows the results for the system with q = 3 pm , g = 4.5 pm and r2 = 7.5 pm . In this case the best fit parameters were: c0 = 8 x10_1° mol cm_3 , k0 = 0.25cms_1, h = 117 s, t2 = 171s , a = 0.2 , P = 0.12 . Fitting results for the
21
third experimental situation are depicted in Fig. 8. The geometry of the system in this case is represented by the following parameters: q = 4 pm, g = 3.5 pm and Г2 = 7.5pm . The best fit parameters are: Co = 1.7 x10_9molcm_3 , k0 = 0.2cms_1, tj = 205 s , t2 = 278s , a = 0.045 , P = 0.6 .
Fig. 7. Fitting results for the system with the dimensions: q = 3 pm, g = 4.5 pm and Г2 = 7.5 pm
Fig. 8. Fitting results for the system with the dimensions: q = 4pm , g = 3.5pm and Г2 = 7.5pm
As can be seen the fitted curves obtained with and without taking into account the homogeneous kinetics are almost indistinguishable. Such behaviour is as expected since the concentrations of NO at the neuron surface are very low (of the order of a few mM or less) and the rate constant of the homogeneous interaction is not sufficiently high to appreciably alter the concentration field of NO°. Thus this means that the influence of the considered homogeneous process will affect the response only at comparatively long times. This was observed in our numerical results. The relative difference in the computed currents obtained with and without the homogeneous term in eq. (7) (see Figs. 68) is below one percent before the current peak. However in the decaying part of the current curve the relative difference increases and, after this wave, reaches values from 5% (the case depicted in Fig. 6) up to 20% (in the third case, Fig. 8).
The neuron activity functions, k(t), with the parameters indicated above for all three cases are shown in Fig. 9.
22
Fig. 9. Neuron activity functions for the three systems shown in Figs. 6-8
5. Conclusions
In this paper, we have developed mathematical models describing the process of nitric oxide release by a stellate neuron and solved them using the ADI method in the conformally transformed space. The developed models represent the diffusion-reaction process which takes place in a slice of the rat’s brain when a patch-clamped neuron is activated using a micropipette. The proposed method of numerical simulation permits accurate and efficient investigation of this process, which in turn allows the solution of the inverse problem of determination of the neuron activity function parameters and the maximum concentration of nitric oxide at the surface of active neuron.
Acknowledgement. We thank Dr. Manon Guille for providing experimental data used in this paper. Dr. Oleinick thanks Maine de la ville de Paris for the post-doctoral fellowship. In Paris, this work was supporting in part by CNRS, ENS (BQR), UPMC and the French Ministry of Research (UMR 8640). In Kharkov, this work is a part of project 158-1.
References: 1. C.S. Roy, C.S. Sherrington // J. Physiology (London) V.11, 1890, pp. 85-108. 2. KCaesar, L.Gold, M. Lauritzen // Proc. Natl. Acad. Sci. USA, V. 100, 2003, pp. 42394244. 3. К. Thomsen, N. Offenhauser, M. Lauritzen // J. Physiology (London) V. 560, 2004, pp.181-189. 4. C. Amatore, S. Arbault, Y. Bouret, B. Cauli, M. Guille, A. Rancillac, J. Rossier // ChemPhysChem, in press. 5. V.L. Pogrebnaya, A.P. Usov, A V. Baranov, A.I. Nesterenko, P.I. Bezyazychnyi // J. App. Chem. USSR. V. 48(5), 1975, pp. 1004-1007. 6. P. C. Ford, DA. Wink, D.M. Stanbury // Febs Lett. V. 326, 1993, pp.1-3. 7. C. Amatore, A.I. Oleinick, I.B. Svir// J Electroanalytical Chemistry. V. 553, 2003, pp. 49-61. 8. C. Amatore, A Oleinick, I. Svir// J Electroanalytical Chemistry. V.564, 2004, pp. 245-260. 9. C. Amatore, A Oleinick, I. Svir// Electrochemistry Communications. V. 5, 2003, pp. 989-994. 10. A Oleinick, C. Amatore, I. Svir // Electrochemistry Communications. V. 6, 2004, pp. 588-594.
Поступила в редколлегию 03.09.2005
Рецензент: д-р техн. наук, проф. Хаханов В.И.
Oleinick Alexander Igorevich, Cand. Tech. Sci., Senior Scientist of the Mathematical and Computer Modelling Laboratory KNURE. Scientific interests: mathematical modelling and numerical methods in mathematical physics, programming. Address: KNURE, Kharkov, 61166, 14 LeninAvenue. Presently a post-doctoral researcher at Ecole Normale Superieure, 24 rue Lhomond, 75231 Paris Cedex 05, France.
Amatore Christian Andre, Professor of Ecole Normale Superieure (ENS), Academician of French Academy of Science, Director of Chimie at ENS. Scientific interests: applied mathematics, physics, analytical chemistry, biochemistry, physical chemistry. Address: ENS, 24 rue Lhomond, 75231 Paris Cedex 05, France.
Svir Irina Borisovna, Chief Scientist, Dr. Tech. Sci., Head of the Mathematical and Computer Modelling Laboratory KNURE. Scientific interests: mathematical modelling and numerical simulation in science. Address: KNURE, Kharkov, 61166, 14 Lenin Avenue. Email: <svir@kture.kharkov.ua>.
РИ, 2005, № 3
