Analysis and simulation of some ODE and PDE models used in neuroscience

Facultad de Matemáticas Departamento de Ecuaciones Diferenciales y Análisis Numérico Trabajo de Fin de Grado Analysis and simulation of some ODE and PDE models used in neuroscience Autor: Alejandro Barea Moreno Dirigido por: José Antonio Langa Rosado Antonio Suárez Fernández Curso 2022/2023 A mamá, a papá y a mi siempre presente amigo Javi. Brain,n. An apparatus with which we think that we think. Ambrose Bierce, The Devil’s Dictionary Alejandro Barea Moreno ODE and PDE models used in neuroscience Index Part 1. The neuron: ion diffusion and action potentials 6 Chapter 1: The Hodgkin-Huxley equations 6 1.1. Anatomyofaneuron................................. 6 1.2. Ion Diffusion and Nernst potential . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3. Theequivalentcircuit................................. 9 1.4. Hodgkin-Huxley equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.4.1. The Problem with the equivalent circuit . . . . . . . . . . . . . . . . . . 12 1.4.2. Voltage-gated channels . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.4.3. Hodgkin-Huxley equations and action potentials . . . . . . . . . . . . . . 15 1.4.4. Quick detour: Hopf bifurcations . . . . . . . . . . . . . . . . . . . . . . . 18 1.4.5. Fixed Points, Limit cycles and Bifurcations . . . . . . . . . . . . . . . . . 19 1.4.6. Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Chapter 2: Reduction to a 2D model 25 2.1. Simplifications of the Hodgkin-Huxley model . . . . . . . . . . . . . . . . . . . . 25 2.1.1. Simulations .................................. 25 2.2. The INa,p `IKmodel................................. 27 2.3. ClassIandClassII.................................. 32 Chapter 3: 1D models 33 3.1. Further simplification of the Hodgkin-Huxley system . . . . . . . . . . . . . . . 33 3.2. LeakyIntegrate-and-fire ............................... 33 3.3. Quadratic Integrate-and-Fire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4. θ-model ........................................ 40 Part 2. When neurons become noisy: the Fokker-Planck approach. 42 Chapter 4: Stochastic Calculus 42 4.1. Randomrealvariables ................................ 42 4.2. Basics of stochastic processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.3. Markovprocesses ................................... 45 4.3.1. Wienerprocesses ............................... 45 4.4. Langevin, Itô, Fokker and Planck . . . . . . . . . . . . . . . . . . . . . . . . . . 46 Chapter 5: Noisy neurons 49 5.1. Populationactivity .................................. 49 5.2. Mean-Fieldargument................................. 50 3 5.3. Boundary conditions, uniqueness and existence . . . . . . . . . . . . . . . . . . . 52 5.3.1. Leaky Integrate-and-fire . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.3.2. θ-model .................................... 54 5.4. Analysisofthesolutions ............................... 55 5.5. StationarySolutions ................................. 58 5.5.1. Leaky Integrate-and-fire . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.5.2. θ-model .................................... 60 5.6. Simulations ...................................... 60 5.6.1. Leaky Integrate-and-fire . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.6.2. θ-model .................................... 62 Conclusions 64 Appendix 65 References 66 Resumen El objetivo de este trabajo es introducir al lector en algunos modelos de EDOs y EDPs que permiten estudiar el comportamiento de la dinámica en el cerebro. Se partirá de algunas hipótesis físicas del comportamiento de las neuronas para deducir modelos que pasarán a ser estudiados cualitativamente y resueltos numéricamente. Posteriormente, se utilizarán aproximaciones para reducir la dimensionalidad y mejorar el análisis para, finalmente, conseguir entender el comportamiento de redes neuronales bajo estímulos estocásticos mediante el uso de la ecuación de Fokker-Planck y de la teoría de ecuaciones diferenciales estocásticas. Abstract The aim of this work is to introduce the reader to some models of ODEs and PDEs that allow us to study the behaviour of the dynamics in the brain. We will start from some physical hypotheses of the behaviour of neurons to deduce models that will be studied qualitatively and solved numerically. Subsequently, approximations will be used to reduce dimensionality and improve the analysis in order to finally understand the behaviour of neuronal networks under stochastic stimuli through the use of the Fokker-Planck equation and the theory of stochastic differential equations. Alejandro Barea Moreno ODE and PDE models used in neuroscience Part 1. The neuron: ion diffusion and action potentials We will dedicate this first part to understand the voltage of neurons. In Chapter 1, we will consider the physical properties of the cellular medium in order to obtain a set of equations, which will then be used to produce simulations. We will simplify these equations in Chapter 2 to a 2D system which is analytically tractable. Finally, in Chapter 3, the 2D system will be further simplified to obtain the Integrate-and-Fire models, single equations which will prove to be very useful in the second part of this work in the end. To understand the work it is necessary to have previous knowledge of the theory of differential equations and their qualitative analysis. The theorems and definitions that will be used throughout the work and that are assumed to have been studied along the degree in Mathematics will be included in an appendix at the end. Chapter 1: The Hodgkin-Huxley equations Alan Hodgkin and Andrew Huxley, two British researchers at Cambridge University, were awarded the Nobel Prize of Medicine in 1963 due to their articles published in 1952 regarding the propagation of action potentials in a squid’s axon. They developed a mathematical framework to understand the change in the potential between the interior and the exterior of cells that have since opened the field of Computational neuroscience. In this first Chapter, we will follow the reasoning they used to write their articles. 1.1. Anatomy of a neuron Fig. 1: Parts of a neuron [1]. 6 Alejandro Barea Moreno ODE and PDE models used in neuroscience First, it is necessary to know some of the biology underlying the processes that will be analyzed. Starting from the basics, a neuron is the basic cell of the nervous system, responsible for the homeostasis, which is the equilibrium and coordination of the different parts of the body. As the figure 1 shows, the neuron is not a homogeneous cell as it has very differentiated parts. Neurons communicate using action potentials which are waves of voltage that travel across the membrane to the synaptic terminals. When they arrive there, the synaptic vesicles release neurotransmitters, chemicals that arrive at other neurons’ dendritic tree and generate post-synaptic potentials which are new action potentials in the activated cells. 1.2. Ion Diffusion and Nernst potential In the cellular medium, there are a lot of ions of different chemical species: Na`, K`, Cl´,etc... We can quantify the difference in voltage generated by the imbalance in the number of ions inside and outside the cell: VM“Vin ´Vout when this voltage surpasses a threshold, an action potential is started. The underlying mechanism of this phenomenon will be thoroughly explained throughout this chapter. Ions move across the cell’s membrane through the so-called ion channels which are proteins that allow the passing of some selective ions. There are hundreds of channels, but we will focus on the sodium and potassium ones, as they are the most important to determine the cell’s voltage. Fig. 2: Ion movement through a cell’s membrane [2]. These ions are subjected to two different forces which will determine their movement through the membrane: diffusion and electric drift. The density of current due to the action of the electric field on the ions (drift) is easily obtained using Ohm’s Law:  jdrift “qt E“ ´c µq∇V(1) 7 Alejandro Barea Moreno ODE and PDE models used in neuroscience where qtis the total charge of the ions travelling through the channel (C), cis the concentration (kg ¨m´3), µqthe electrical mobility (m2¨V´1¨s´1). Diffusion is a phenomenon originating from the Brownian (random) motion of the particles in the liquid. When observing a gradient of concentration in a chemical species, the spontaneous motion of the particles macroscopically is towards the zone of less concentration. The action of special enzymes called bombs will maintain this concentration gradient, but this will not interfere with our analysis. Mathematically this is modelled by the first of Fick’s Laws:  jdiff “ ´D∇c(2) where cis the concentration and Dis the diffusion constant (m2¨s´1). Einstein found that this constant is related to the mobility of the ions (µq): D“µqkBT q“µqkBT ze (3) where kB“1.38 ¨10´23 m2¨kg ¨s´2¨K´1is the Boltzmann’s constant, Tthe temperature (K) and e“1.6¨10´19 Cthe elementary charge, which is the net charge of a proton or an electron. The zin this equation is the valence of the ions, which is the number of net charges that it contains, as the ion can have an excess or a defect of electrons, zis an entire number. Now, if we assume there is no interaction between ions in a channel and that this channel is approximately unidimensional, we can sum (1) and (2) and use (3) to obtain: jT“ ´µqcdV dx ´µqkBT ze dc dx.(4) Which is the total density of flux of ions through the channel. As we have assumed that the channel is unidimensional, the gradients are reduced to derivatives in the spatial variable xwhich represents the length of the channel. It will be interesting to know the value of the voltage needed to stop the flow of a particular ion, which will be easily obtained by equating (4) to zero and then integrating between the inside and outside of the cell: Ei“kBT zieln ˜cin,i cout,i ¸. This is the Nernst Potential of the ion i. This potential represents the membrane voltage at which the flux of that particular ion is stopped, which will be illustrated in the next section. 8 Alejandro Barea Moreno ODE and PDE models used in neuroscience Fig. 6: Evolution of the functions involved in equations (14). a) Steady-state of the gating variables. b) Time constants of the gating variables. The input must be expressed in (mV ) and the output is in (ms´1). These functions are all continuous and non-negative. This is easily seen for some, but for others, simply applying L’Hôpital’s rule allows us to define them at their conflicting points so that they are continuous. With this Hodgkin and Huxley modelled the variation of the conductances: gK“gKn4, gNa “gNa m3h. (15) The exponents have biological meaning, the exponent 4 in the potassium is due to the four particles needed to be attached to the membrane to open the gate. The hvariable in the sodium gate is an inactivation variable (biologically a blocking particle) and the exponent 3 over the mare the 3 particles needed for the gate to open [4]. 1.4.3. Hodgkin-Huxley equations and action potentials With these gating variables, the circuit equation (10), the equations (14) for m,hand n and (15) we have finally deduced the Hodgkin-Huxley equations: CM dVM dt “ ´gKn4pVM´EKq´gNa h m3pVM´ENaq´gLpVM´ELq`Iptq(16) dm dt “m8pVMq´m τmpVMq(17) dn dt “n8pVMq´n τnpVMq(18) dh dt “h8pVMq´h τhpVMq.(19) 15 Alejandro Barea Moreno ODE and PDE models used in neuroscience These equations are only valid when, by using techniques which are not important as of now, we attain a homogeneous voltage throughout the cell. This system of coupled non-linear equations is non-analytically solvable. However, these functions are continuous and have continuous partial derivatives so we can assure that for each initial condition, they have a solution and this solution is unique, we know this because of theorems 0.6.1 and 0.6.2. Lemma 1.4.2. Given an initial condition ϕpt0q “ pV0, m0, n0, h0qwith 0ďm0, n0, h0ď1 and the steady-states and time constants of the gating variables are bounded, then the equation (13) has a unique solution ϕptq “ pϕVptq, ϕmptq, ϕnptq, ϕhptqq that exists in rt0,8q and 0ď ϕmptq, ϕnptq, ϕhptq ď 1,@tP rt0,8q. Proof: We must prove now that the solution is defined for each tět0“0. For this, we must find a bound the positive trajectory of the system. We know that the solution is bounded in the m,nand hdimensions (Lemma 1.4.1). Let ϕptq“pϕVptq, ϕmptq, ϕnptq, ϕhptqq the solution of (16) given the initial condition ϕpt0q “ pV0, m0, n0, h0qthen by using the integral formulation of the Cauchy Problem we have: ϕVptq “ V0`żt t0p´gtotpsqpϕVpsq´Etotpsqq`Ipsqqds. (20) Where: gtotptq “ gNaϕhptqϕmptq3`gKϕnptq4`gL.(21) Etotptq “ ENagNaϕhptqϕmptq3`EKgKϕnptq4`ELgL gtotptq.(22) By using what we proved in Lemma 1.4.2 and the fact that gtotptq ą gL,@t, we know that (21) and (22) are both bounded. Now we subtract Etotptqfrom both sides and apply absolute values to both sides of the equation, obtaining, after some algebra: |ϕVptq´Etotptq|ď|V0|`|Etotptq|`∥I∥8pt´t0q`żt t0 gtot|ϕVpsq´Etotpsq|ds. (23) Now, let us consider that our solution only stays finite in [t0,t˚) and therefore in [t0,t1] with t1“t˚´ϵand that IPL8pRq, then we can do: |V0|`|Etotptq|`∥I∥8pt´t0q ď M. (24) 16 Alejandro Barea Moreno ODE and PDE models used in neuroscience By applying (24) in (23) we can use Gronwall’s Lemma (0.6.1) to obtain: |ϕVptq´Etotptq|ďMeşt t0gtotpsqds ďMept˚´t0q∥gtot∥8.(25) As gtotptqand Etotptqare bounded, for each finite time t˚, the trajectory remains bounded and therefore can be extended. This concludes our proof that our solution exists in rt0,8q. ˝ Looking at the equations and graphs we have already shown above with the steady states and time constants (which are not constants, but we will keep this nomenclature), we can qualitatively analyse what is going to happen before we go into more detail. The stimulus that will allow the neuron’s potential to overcome the threshold will be given by the intensity represented either by a current injected by the experimenter or by the action potentials of other neurons arriving via the synapse. If the threshold is not reached, there will be a small depolarisation followed by a re-establishment of the equilibrium potential, corresponding to the fixed point of the unexcited system. However, if the value that allows the potential to fire is reached, what will happen is that the variable mwill rapidly go to steady state, due to the fact that τmis lower than the other time constants, causing the sodium term to be the most important at first and the membrane to depolarise (increase its potential) to values close to the sodium Nernst potential, reducing the term pVM´ENaq. By the time this happens, the variable hdecreases in value as nincreases, giving more importance to the potassium term and causing it to hyperpolarise (decrease in potential) to values close to the potassium Nernst potential. Once this occurs, ndecreases and the equilibrium state is reached again, with no more action potentials being released until a new stimulus is imposed. This is summarised in the following picture: Fig. 7: Action potential explained by means of the dynamical Hodgkin-Huxley system [5]. 17 Alejandro Barea Moreno ODE and PDE models used in neuroscience 1.4.4. Quick detour: Hopf bifurcations We will study the behaviour of the dynamical system of Hodgkin and Huxley by using elements from the qualitative theory of differential equations, mainly Theorem 0.6.3. However, there is something, which is not studied in preliminary courses of differential equations which must be discussed before going further and that is theory of structural stability. If we have a dynamical system dependent on a parameter, as is our case here with I(we will ignore the variability with respect to Nernst conductances and voltages), it is possible that the behaviour of some equilibria or limit cycles of the system changes according to the value of this parameter, i.e. the system does not remain structurally stable. This is what is behind action potentials, which are periodic phenomena generated by the approach of a trajectory to a closed orbit that will act as a limit cycle. These potentials do not occur for any value of the parameter I, but for those in which the conditions are conducive to what we have just explained. First of all, I would like to mention that the phenomena we will be studying will arise from the change in the local behaviour of the fixed points of the system and that we can infer this behaviour using the linear version of our dynamical system in the vicinity of these points thanks to the following Theorem [6]: Theorem 1.4.1 (Hartman-Grobman Theorem).In a vicinity of a hyperbolic fixed point, a dynamical system and its linearized version are topologically equivalent. The definition of hyperbolic fixed points is stated in the appendix (0.6.1). Without entering into much detail about what it means to be topologically equivalent, we will say that two dynamical systems are topologically equivalent when you can find a continuous function that maps trajectories from one system to the other. The bifurcation with which we will be working in this chapter will be the Hopf bifurcation which occurs when two complex eigenvalues of the Jacobian matrix of the system cross the imaginary axis and a fixed point changes from being stable to unstable or vice versa. Hopf bifurcations come in two flavours: subcritical and supercritical. In subcritical bifurcations, there is a point before the bifurcation when an unstable limit cycle and a stable one, both of them around a stable equilibrium, appear. Therefore, trajectories starting between the two cycles will become periodical and the ones starting inside the unstable limit cycle will eventually end at the fixed point. When the threshold is surpassed the unstable limit cycle and the fixed point fuse, making the latter unstable and generating oscillating orbits. These systems present hysteresis, because after lowering the parameter below the Hopf bifurcation, oscillations may not stop until achieving the point we talked about before where the two cycles appear. 18 Alejandro Barea Moreno ODE and PDE models used in neuroscience Fig. 8: Behaviour of a subcritical Hopf bifurcation generated by the parameter µ. The unstable limit cycle is represented with the dashed line which we see disappearing and turning the point unstable (white) as µą0[3]. In supercritical bifurcations, the oscillations occur little by little after the threshold is surpassed. The stability is lost as from the point appears a stable limit cycle. Fig. 9: Behaviour of a supercritical Hopf bifurcation generated by the parameter µ. A stable limit cycle appears as µą0[3]. 1.4.5. Fixed Points, Limit cycles and Bifurcations We will consider Iptqas a constant input and study the fixed points depending on its value. By equating Hodgkin-Huxley equations to zero, we obtain that the fixed points will be the solutions to the equation: I“gKn8pVMq4pVM´EKq`gNah8pVMqm8pVMq3pVM´ENaq`gLpVM´ELq “ FpVMq.(26) So the problem is now simplified, as it only consists in analyzing the number of zeros that the function Fhas. We have omitted the CMbecause in the units we are using its value as 1. We will be using the data from table 1 for our simulations, as the derivative of (26) will be difficult to discuss analytically, I have done a numerical experiment with MATLAB which shows that the function is monotone. It is also easy to see that the function is surjective and therefore we have a unique solution for each I, which we will call V˚ MpIqfrom now on. We now must analyze the local stability of the system by linearizing it and studying the sign of its eigenvalues. The jacobian matrix of the system of Hodgkin and Huxley has the following 19 Alejandro Barea Moreno ODE and PDE models used in neuroscience form: JpV˚ MpIqq “ ¨ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˝ ´BF BV´BF Bm´BF Bn´BF Bh m1 8 τm 1 τm 0 0 n1 8 τn 01 τn 0 h1 8 τh 0 0 1 τh ˛ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‚ (27) the analysis of the sign of the eigenvalues of the matrix (27) can easily be done numerically (which we will soon be doing) but an in-depth analysis requires a lot of mathematical tools which we will not show. However in [7], Troy shows that there exist two Hopf bifurcations in the system depending on I. We can show this numerically in a graph: Fig. 10: Evolution of the local stability of the fixed point of Hodgkin-Huxley equations. We represent the value V˚ MpIqas a function of I. Red points are unstable and black points are stable. Our estimations locate the two bifurcations at: I1 Hopf “8.44µA ¨cm´2, I2 Hopf «163.37µA ¨cm´2 1.4.6. Numerical Simulations Now, with the simulations. We start at equilibrium with I“0. (Vp0q“´64.98mV, mp0q “ 0.05, np0q “ 0.32, hp0q “ 0.60). First of all, we will use short current pulses which imitate the behaviour of a Dirac delta function: 20 Alejandro Barea Moreno ODE and PDE models used in neuroscience Fig. 11: Response of a neuron predicted by the H-H equations under a Dirac delta-like input fired at 20 ms. As we see, the action potential is an ’all of nothing’ phenomenon, as we had explained before. The threshold we have calculated before is not valid in this type of pulse because it is very sensitive to their duration. We can understand these Dirac delta pulses as finite and instantaneous jumps in coordinates in the 4D phase plane. When the jump is big enough, the flow of the system starts an action potential which then finishes in the stable equilibrium we calculated before. Now, if we consider a constant current with greater amplitude than the one we found as the trigger in the latter section, we obtain the following behaviour of the different variables: Fig. 12: Response of a neuron predicted by the H-H equations under a step current injected and stopping at the times marked by the red lines. I“10µA ¨cm´2. 21 Alejandro Barea Moreno ODE and PDE models used in neuroscience We observe repetitive firing, something which is known as bursting in the world of neuroscience. The behaviour that we expected is what we actually see, the first gating variable to fire is the mone whose peaks are narrower than the others because of the smaller magnitude of its gating variable. Also, I have marked with green dashed lines the Nernst potentials of sodium and potassium and we see how it approaches both but it comes nearer the potassium one („ -77 mV ). When the step current is stopped, the system swiftly approaches rest. The shape of the peaks of the nand hgating variables is interesting because we see that one is almost the exact opposite of the other, we will make use of this fact in the following chapter to reduce the dimension of the HH system. Another thing that comes to your mind when seeing this picture is that the first action potential is much higher than the rest, this is due to the solution approaching the limit cycle in spiral motion. This can be easily observed in figure 13: Fig. 13: Phase-plane portrait of the former situation. The red point is the start of the orbit. The black arrow signals the first action potential as it approaches the limit cycle and the red arrow the return to the initial state of rest after firing 7 times (doing 7 turns). This is proof of the usefulness of this model because this is something that is actually observed in neurons when they fire several action potentials, which is the bursting we mentioned before. Now, we can calculate the frequency of these action potentials as a function of I, as it is characteristic of the type of model that we have, as we will see later, which gives us the following graph: 22 Alejandro Barea Moreno ODE and PDE models used in neuroscience Fig. 14: f´Igraph for a H-H neuron. Due to the computational difficulty of calculating these graphs, we cut before the second bifurcation occurs. The important detail of this one is the discontinuity that occurs near the first bifurcation, which we will see is characteristic of Hopf-type bifurcations. The fact that oscillations begin before the value of the bifurcation is a hint that maybe we have a subcritical Hopf bifurcation. We can do a little numerical experiment to see another hint that the first bifurcation may be subcritical. We will model the intensity as a decreasing ramp and this is what we observe: Fig. 15: Response of the HH equations under an input Iptq “ 20 ´0.1¨t µA ¨cm´2. The red dashed line indicates the point at which Iptqequals the bifurcation value. 23 Alejandro Barea Moreno ODE and PDE models used in neuroscience The hysteric behaviour is typical of subcritical Hopf bifurcations. However, this is by no means a rigorous proof, even though the literature tells us that it is in fact, subcritical [2], [7]. For the second bifurcation, we can show that the amplitude of the oscillations begins to decline long before it is surpassed: Fig. 16: Peak to peak voltage of the signals as a function of the applied intensity. Therefore, it seems that this second bifurcation is supercritical. All in all, we have successfully obtained the Hodgkin-Huxley equations, the most fundamental differential system in neuroscience. It is worth saying that we can extend this system to include the influence of each one of the channels that are present in a cell’s membrane. However, this will increase notably the computational effort to calculate the solutions. 24 Alejandro Barea Moreno ODE and PDE models used in neuroscience Fig. 24: Limit cycle that generates the periodic firing in the INa,p ´IKmodel with I = 6 µA ¨cm´2. talked about before when we guessed that we are experiencing a SNIC bifurcation. We can see the f´Idiagram of these neurons and we find something interesting: Fig. 25: f´Igraph for a INa,p ´IKneuron. Results have been recorded over a 1s pulse. In this case, we do not observe a discontinuity, this is a consequence of the internal workings of a saddle-node bifurcation, where the frequency dependency is proportional to ?I´ISNIC. In the next chapter, we will work with another model that shows the same behaviour and obtain this relation for every model that has the same type of bifurcation. The fact that the frequency starts from 0 and is continuous makes it possible to adjust the distance between subsequent action potentials. 31 Alejandro Barea Moreno ODE and PDE models used in neuroscience 2.3. Class I and Class II We have seen two types of behaviours in the f´Icharacteristic of our models, these were categorised by Hodgkin in two different classes: Class I: They can start firing at arbitrarily low frequencies and the time between peaks is short. Class II: They have a discontinuity in their f´Igraph and the recovery time between peaks is longer. Neurons that do not belong to any of these two groups are viewed as abnormal. The remarkable thing is that the shape of these curves depends on the type of bifurcation that serves as a threshold for their action potentials: models with Hopf bifurcations generate class II graphs while models with Saddle-node bifurcations generate class I f´Icharacteristics. Neurons of different parts of the brain behave in one way or another, this classification is exhaustive as neurons that do not correspond to either of the categories are considered anomalous. 32 Alejandro Barea Moreno ODE and PDE models used in neuroscience Chapter 3: 1D models In this chapter, we will be studying 1-dimensional models specifically designed for their application in neural networks. We will start from the famous Leaky Integrate-and-fire model and conclude with the ring models used for neurons in the visual cortex. 3.1. Further simplification of the Hodgkin-Huxley system In the former chapter, we studied two models whose orbits follow the dynamics given by systems of equations of the form: τM dVM dt “FpVM, wq`RMIptq τw dw dt “GpVM, wq where τM“RMCM. In the actual experience, it is usually seen that τwąą τMand thus the derivatives in the phase plane are almost horizontal. Therefore it is entirely justified to suppress the second equation and obtain the following: τM dVM dt “FpVMq`RMIptq.(29) This is the general equation of Integrate-and-fire models. It is obvious that a local solution given an initial condition Vpt0q “ V0exists and is unique if Fis sufficiently regular, which will be the case in the following cases. Before starting the analysis of the models, from now on we will not give too much importance to the units of the magnitudes because, in this type of model, abstractions are used that go beyond the concrete physics of the models. This will not be a big problem since units have not played an important role up to now. 3.2. Leaky Integrate-and-fire The leaky integrate-and-fire is a very simple model in which the F function in equation (29) is a linear one: τM dVM dt “ ´pVM´Vrestq`RMIptq.(30) This equation is obtained through a circuit like the one in (10) but with only one channel with resistance RMand then τ“RMCM. This model is called leaky because we see that the only current that we consider is the "leak" one, the one with constant conductance. We can 33 Alejandro Barea Moreno ODE and PDE models used in neuroscience easily find a general solution given an initial condition Vpt0q “ V0: VMptq “ V0exp ˆ´t´t0 τM˙`Vrest «1´exp ˆ´t´t0 τM˙ff` exp ˆ´t τM˙ CMżt t0 Ipsqexp ˆs τM˙ds. (31) This system has the same response that we saw in the first chapter when we approximated the membrane by an equivalent circuit. The interesting thing that these models do and why they are useful is that they serve us to sum (integrate) the inputs that other neurons send when they fire an action potential, modulated by the Green function of the system which is, in this case, exp ˆs τ˙. The fact is that, as we had seen with the equivalent circuit in Chapter 1, this is not the behaviour we want. With this model let alone like this we will not observe repetitive firing. It only has a fixed point which is always attractive as the Jacobian of the system is always ´1, therefore all of the trajectories will be attracted to it. The solution to this problem is to insert a manual threshold. If the membrane potential arrives at a value Vth, then it is reset manually in Vreset ăVth after waiting a time which we call the refractary time tref , analogue to the one we saw in actual action potentials. For the sake of simplicity and mathematical cleanliness, we will only consider the cases such that Vreset “VRand tref “0. This gives us the following properties for the model: Property 3.2.1. When considering a constant current input Iptq “ I, the system has a current threshold for firing given by the expression: Imin “Vth ´VR RM . Proof: Let us consider a constant current to determine the threshold. The steady state of the system is given by: 0“ ´pVp8q´VRq`RMI. (32) For a system to fire, the membrane potential must arrive at Vth, the most extreme case is that the system arrives at this threshold at infinity, therefore, substituting Vp8q with Vth in (32) we obtain the result. ˝ Property 3.2.2. If V0“VR, the frequency of firing is given by the equation f“1 τMln ˆ1`Imin I´Imin ˙,(33) 34 Alejandro Barea Moreno ODE and PDE models used in neuroscience which behaves linearly at high currents. Proof: If the current is constant and VMp0q “ V0“VR, the equation (31) turns into Vptq “ VR`RMIp1´e´t τMq, because tref “0and Vreset “VR, we know that the period of the signal is the time that it takes for it to go from the reset to the threshold and can be thus calculated: Vth “VR`RMIp1´e´T τMq where, doing some algebra and using Property 1, we obtain: f“1 T“1 τMln ˆ1`Imin I´Imin ˙. Then, for Iąą 1, we can use the Taylor expansion of the logarithm: lnp1`xq “ x`opxq. And we obtain: fIąą1«I´Imin τMImin which is linear. ˝ Corollary 3.2.1. The LIF model is class I. Using the following data: RM0.125 kΩ τM0.125 ms VR-65 mV Vth 40 mV Table 3: Data for the simulations with the Leaky Integrate-and-fire model. Using the formula given by equation (3.2.1), we can estimate the value of the minimum current with these parameters: Imin “840µA. We can do a simple simulation to show the behaviour of the model: 35 Alejandro Barea Moreno ODE and PDE models used in neuroscience Fig. 26: Simulation of the Leaky Integrate-and-fire model under a constant current stimulus. In orange if IăImin and in blue IąImin. The red dashed line indicates the position of the threshold. When the current is not smaller than the threshold (Imin “840µAqthe system does not fire, as predicted, but when it is greater, the repetitive firing cannot be stopped. In terms of the phase space (or in this case phase-line), it is easy to see why this happens. The steady-state if IąImin will be over Vth and, as we said before, will be attractive. When the points approach the equilibrium but surpass the threshold, the system sends the trajectory back to the starting point.1 3.3. Quadratic Integrate-and-Fire The function FpVMqwhich appears in (29) has a non-linear behaviour which cannot be understood solely by the leaky integrate-and-fire. This is what motivates the introduction of the Quadratic Integrate-and-Fire: dVM dt “apIptq´I1q`bpVM´V1q2(34) where ahas units of (F´1) and b(V´1¨s´1) and with I1, V1are all positive parameters. In the case in which Iptq “ I, which is the one we will be working from now on, this equation can be analytically solved, but before finding the solution, let us study something about the behaviour of fixed points: 1Like Sisyphus, it must repeat the process over and over again without reaching the top of the mountain. 36 Alejandro Barea Moreno ODE and PDE models used in neuroscience Property 3.3.1. The equation (34) has two fixed points if I1ąIand none in the complementary. Their coordinates will be: V˚ M“V1˘ca bpI1´Iq.(35) Proof: Using equation (34), and developing the equation for fixed points we obtain: pVM´V1q2“a bpI1´Iq. Because aand bare positive, we know that this equation only has real solutions if I1ąI, and the coordinates of such fixed points will be: VM“V1˘ca bpI1´Iq. ˝ Property 3.3.2. If I1ąI, one of the fixed points is stable and attractive and the other is unstable. Proof: The derivative of the system with respect to VMis: JfpVMq “ 2bpVM´V1q.(36) Substituting (35) into (36) we obtain: JfpV˚ Mq“˘2bca bpI1´Iq. Therefore the one to the left will be attractive and stable and the one to the right unstable. ˝ Corollary 3.3.1. The system has a saddle-node bifurcation in I“I1, VM“V1. Property 3.3.3. When I1ăIor I1ąIand VMpt0q ą V1`aa bpI1´IqV is monotonically increasing. Lemma 3.3.1. When I1ăIor I1ąIand VMpt0q ą V1`aa bpI1´Iqthe solution explodes in finite time. Otherwise, it stays finite. Proof: The proof of the second affirmation, that is, the case in which I1ąIand Vpt0q ă V1`aa bpI1´Iq, is trivial seeing what we just studied about fixed points. For the remainder of the proof, t0“0. 37 Alejandro Barea Moreno ODE and PDE models used in neuroscience Equation (34) is a Ricatti differential equation, which we can reduce to: d2x dt2“γ2x, γ “a´abpI´I1q.(37) With the change of variables VM“V1´1 xb dx dt , the solution to equation (37) is: xptq “ Aeγt `Be´γt. Therefore: VMptq “ V1´γ b˜A Beγt ´e´γt A Beγt `e´γt ¸. Imposing the initial condition VMp0qwe obtain: A B“1´D 1`D, D “bVMp0q´V1 γ. For the solution to explode in finite time, there must be a value of t, let us say Twhere: 1´D 1`DeγT `e´γT “0. Now, we must distinguish between cases. If IąI1, then γis pure imaginary and the former equation can be transformed into: tanpiγTq “ iγ bpVMp0q´V1q which has a solution because the tangent is surjective. If IăI1,γis real and: e´2γT “D´1 1`D. As the exponential is positive the equation has a solution if and only if Dě1which means: Dě1ñVMp0q ě V1`c´a bpI´I1q. Which proves the result. ˝ So, what can we do when our model goes to infinity on a common basis? We consider the times at which the solution does this as firing instants and the reset potential is ´8. Then: Property 3.3.4. The frequency of firing if IąI1fulfils the following law: 38 Alejandro Barea Moreno ODE and PDE models used in neuroscience f9aI´I1. Proof: We know that if IąI1, we can use the equation: tanpiγTq “ iγ bpVMp0q´V1q to compute the period. If we impose VMp0q “ ´8 then: tanpiγTq “ 0ñT“π iγ “π aabpI´I1q. Taking the inverse gives us the desired result. ˝ Using the parameters: a1F´1 b1mV ´1¨ms´1 V10 mV I12 A Table 4: Data for the simulations with the Quadratic Integrate-and-fire model. We can do some simulations to show what we just proved (of course we stop before going to 8, because we cannot treat that numerically): Fig. 27: Simulation of the QIF model with I“2.5A, V p0q “ ´8, which is over the threshold. The red dashed line indicates the position of the threshold. 39 Alejandro Barea Moreno ODE and PDE models used in neuroscience From this we have measured the distance between peaks obtaining an estimate which is near the value that we can calculate using the last property: Texp “4.42 ms «4.44 ms “π aabpI´I1q which makes sense because of the fact that we stopped the simulation before actually going to infinity, so the simulated period is slightly smaller. If I“1ăI1and the initial condition is between the fixed points we find, as expected: Fig. 28: Simulation of the QIF model with I“2A, V p0q “ ´2mV , which is under the threshold. The red dashed line indicates the position of the threshold. 3.4. θ-model Let the following change of variables in (34): VM“V1`c btan ˆθ 2˙ where c is a positive constant (s´1). Then, it can be found that the equation for θptqis: dθ dt “cp1´cospθqq` ab cp1`cospθqqpI´I1q.(38) This is the so-called θ-model. 40 Alejandro Barea Moreno ODE and PDE models used in neuroscience Definition 4.4.1 (Itô integral).Let fptqa piecewise continuous functions, we define the Itô integral (whenever it converges) as: żt t0 fpsqdBs“lim nÑ8 n ÿ j“0 fpτjq∆Bj where τj“tj. In the case of (45), we have that: żt t0 bpXs, sqdBs“lim nÑ8 n ÿ j“0 bpXtj, tjq∆Bj. The Itô integral has a very interesting property: Theorem 4.4.1 (Itô formula).Let Xtbe the solution of an equation as in (45) with bounded coefficients. Let gpx, tq P C2pr0,8qˆRqand let us define: Yt“gpXt, tq. Then Ytis the solution of: dYt“Bg BtpXt, tqdt `Bg BxpXt, tqdXt`1 2B2g Bx2pXt, tqpdXtq2 where we must consider: dt ¨dt “dt ¨dBt“dBt¨dt “0, dBt¨dBt“dt. From this, we can finally deduce the Fokker-Planck equation. Theorem 4.4.2 (Fokker-PFlanck equation).The probability function P1px, tqcorresponding to the solution of the Langevin equation (45) fulfills: BP1px, tq Bt“ ´BpP1px, tqapx, tqq Bx`1 2B2pP1px, tqbpx, tq2q Bx2 for almost all px, tq P Rˆr0,8q. Proof: Let us consider ha function that is only dependent on xand apply Itô’s formula: dYt“dhpXtq dx apXt, tqdt `dhpXtq dx bpXt, tqdBt`1 2 d2hpXtq dx2b2pXt, tqdt. 47 Alejandro Barea Moreno ODE and PDE models used in neuroscience Taking derivatives with respect to time and taking the expected values from both sides we obtain: B BtżhpxqP1px, tqdx “ż˜dh dxapx, tq` 1 2 d2h dx2b2px, tq¸P1px, tqdx. If we do integration by parts on the second side, we obtain: B BtżhpxqP1px, tqdx “ż«´d dxpapx, tqP1px, tqq` 1 2 d2 dx2pb2px, tqP1px, tqqffhpxqdx. Because the function his arbitrary, we can just, assuming enough regularity of P1move the derivative with respect to time into the integral and obtain the result. ˝ We know that the solution of the Fokker-Planck equation is the probability distribution of a Markov process because we have used this hypothesis to deduce it. However, in [10] we can see a Lemma that shows that if the coefficients of the Langevin equation are independent of time, then the process will be time-homogeneous. 48 Alejandro Barea Moreno ODE and PDE models used in neuroscience Chapter 5: Noisy neurons We have introduced some mathematical tools that helped us to deduce the Fokker-Planck equation in a general case. However, how does this apply to neural networks? What are we expecting to obtain? Which are the boundary conditions that we must impose for our solution to have physical sense? In this chapter, we will answer all of these questions and give a physical context to mathematics. 5.1. Population activity There are many, many neurons in the brain („100 billion [13]), each of which evolves according to its own dynamics. In reality, our calculations so far have been somewhat anecdotal and illustrative, but if we want to make predictions about the behaviour of the brain in such a way that they can have something to do with the macroscopic sphere and with the data obtained in diagnostic tests, we cannot restrict ourselves to one or two neurons, we have to talk about a significantly larger set of neurons: Neuronal Populations. Each neuronal population comprises neurons that have similar properties and respond in the same way to the same stimuli, which will mathematically translate into their dynamics following the same differential equations. An example of a neuronal population is neighbouring cells in the visual cortex that respond favourably to a given orientation. In this case, the main orientation would be the discrimination criterion that would allow us to define the boundaries of the population. Let us imagine now that we have a set of N neurons that are part of the same population. Each of them will have its own voltage and will communicate with each other by sending action potentials. If the dynamics of each cell are governed by a system of k equations and we take into account that they are not independent of each other, in order to know the state of the system, that is, the voltages of each of the neurons, we would need to solve a system of k N equations, and knowing that N is usually a huge number, the reader will be able to deduce that this is practically an impossible task. This situation is very similar to what happened at the dawn of statistical mechanics. When one wants to study the kinetics of gases, one realises that it is impossible to solve Hamilton’s equations for systems with such a large number of particles. Thus, we began to work not with microscopic quantities such as the momentum and position of individual molecules, but with macroscopic quantities such as pressure and temperature, deduced from the probability of the particles occupying specific volumes. We will perform an analogous treatment to this, where the spatial volumes will be replaced by abstract volumes in the dimension of the membrane voltage. Moreover, as we have been introduced in the previous chapter, our treatment will be time-dependent. But what is the macroscopic quantity we want to measure? 49 Alejandro Barea Moreno ODE and PDE models used in neuroscience We want a quantity that gives us information about the functioning of a neuronal population, and a way to characterise it. This will be the firing rate of action potentials in the network or Population Activity. There are several ways of defining this firing rate, which we will denote by the letter A. First, we can think of the most obvious definition: Aptq “ npeakspt, t `Ttestq Ttest (46) where npeakspt, Ttestqis the number of peaks between t and t `Ttest. While this is the most intuitive, it does not really make physical sense, since, in nature, the brain does not wait to estimate the firing rate, because, if it did, phenomena such as reflexes would not exist. But this is where the role of a large number of neurons in our brain comes into play. Indeed, if we average across our population we get a much simpler and quicker result. Thus we get a definitive definition of the magnitude to be studied: Definition 5.1.1 (Population Activity).The activity of a neuronal population with N cells is defined by: Aptq “ lim ∆tÑ0 npeakspt, t `∆tq N∆t.(47) And it is measured in Hertz (Hz). 5.2. Mean-Field argument From now on, we will work with 1D models. We have already mentioned about the humongous number of neurons that inhabit the brain and even though we are treating a limited portion of them, namely a neuronal population, the number is also large enough to make the treatment of a system of N equations, being N the number of neurons in said population, nearly impossible. One of the facts that makes this process very difficult is the inability of differentiating between neurons and thus of differentiating between current inputs. This is why we need to resort to stochastic theory and assume that the input current of a neuron is random: Iptq “ µptq`σptqBt(48) where µptqis the mean, σthe standard deviation and Bta Wiener process or white noise which we defined before. For the mathematics of the following sections to be easier and more satisfactory, we will assume that µptq “ µ, σptq “ σ. Even though a time-dependent analysis would prove interesting, 50 Alejandro Barea Moreno ODE and PDE models used in neuroscience using constant parameters is not a bad approximation based on what is actually observed in experimentation. Now, if we consider a 1D model whose function does not include the time explicitely: dVM dt “FpVM, Iq. If we insert equation (48) into the value of I, we obtain a Langevin equation: dVM“fpVM, µqdt `gpVM, σqdBt. And therefore, using the results proven in the former chapter, we obtain the following Fokker- Planck equation for the Langevin equation (5.2): BPpVM, tq Bt“ ´BpfpVM, µqPpVM, tqq BVM`1 2B2gpVM, σq2PpVM, tq BV2 M .(49) This equation will give us the probability of a neuron’s membrane voltage of having the value VMat an instant of time t. We want to study the ensemble behaviour and for this, we will use arguments inherited from statistical mechanics. We know that the input Iiptqof a neuron is given by equation (48) and this equation does not contain any reference to the index of the neuron, therefore is independent. If we assume that the neurons are identical, as they receive the same input, if we study the stationary scenario, we will deduce that the population activity is equal to the firing rate of an individual neuron. Lemma 5.2.1. If the activity of a neuronal population Aptqis constant, then it is equal to the firing rate of an individual neuron. Proof: Suppose that Aptq “ A0. We know that in a time period T, the N neurons of the neuronal population will produce a number of spikes equal to: nnetwork spikes pT, Nq “ NT A0. If we call the stationary value of the firing rate of a neuron f0, then we know that in the time period T, it will produce: nneuron spikes pTq “ Tf0. Therefore, as all of the neurons are identical, the spatial average that we take in the definition of the population activity (46), we get: 51 Alejandro Barea Moreno ODE and PDE models used in neuroscience A0“f0. ˝ If we now assume that the dynamics of the system are quasi-stationary, that is, at each time equilibrium is reached, we can assume that the activity of the population, which is a magnitude that averages itself, is equal to the firing rate of a neuron, and the Fokker-Planck equation (49) is used to study it. The meaning of this new probability density function is the following: żV2 V1 PpVM, tqdVM“neurons with membrane voltage between V1and V2 N. This application of the so-called mean-field arguments is not very rigorous, but it serves to illustrate the main ideas that lead us to disregard the individual character of the neurons that make up the network. 5.3. Boundary conditions, uniqueness and existence Obviously, we cannot solve this equation as it is, but we need some good boundary conditions. First of all, it is clear that we have to give it an initial profile P0pVMq “ PpVM,0qto work with, and in the next chapter, we will be testing the effects that are generated depending on the type of temporary initial condition we give. On the other hand, we know that what we have in hand is a probability distribution. Let Ω be the spatial domain of the solution (which in our case will be unidimensional), we will have to impose that: żΩ PpVM, tqdVM“1,@tP r0,8s.(50) If this condition is fulfilled by the initial profile, it suffices to prove that the temporal derivative of the integral is zero. It is interesting the choice of writing the Fokker-Planck equation (49) as: BPpVM, tq Bt“ ´BJpVM, tq BVM (51) where JpVM, tq “ fpVM, µqPpVM, tq´ 1 2BgpVM, σq2pVM, tq BVM (52) is a flux of probability. We have dropped the µand σdependency for a neater presentation. 52 Alejandro Barea Moreno ODE and PDE models used in neuroscience This new form helps us understand the fact that the Fokker-Planck equation is indeed a conservation equation, similar to those of charge or fluid mechanics. We will be treating two models stochastically: the Leaky Integrate-and-fire (30) and the θ-model (6), which is why we divide the discussion in two parts. 5.3.1. Leaky Integrate-and-fire Now, we should consider the fact that we will be working with Integrate-and-fire models. In these models, we imposed a threshold for an action potential to be fired. The way to do this was kind of rudimentary: the potential had to be manually reintroduced at the reset voltage Vreset. Here, we can impose the threshold by adding a Dirichlet boundary condition: PpVth, tq “ 0,@tP r0,8q.(53) From now on, the domain of the solution in the voltage dimension will be p´8, Vths. Also, for the probability to be normalisable we need to impose: lim VMÑ´8 PpVM, tq “ 0,@tP r0,8q.(54) This condition, in the context of the Fokker-Planck equation, is called an absorbing boundary condition. This is because all points (neurons) that reach the right-hand boundary, the one defined by Vth, will vanish. Even if P cancels out at this boundary, there is a non-zero flux J and therefore all neurons passing its refractory period will vanish from our probability distribution, to solve this a delta term is introduced to the equation: BPpVM, tq Bt“ ´BJpVM, tq BVM`JpVth, tqδpV´Vresetq.(55) Which means that the probability flux leaving on the far right will be reinserted at Vreset ă Vth which is the reset voltage of our neurons. This term will cause many problems for the numerical solving of the solutions, so we will have to approximate the Dirac delta by a function that resembles it, has integral 1 in the domain of our solution and is smooth enough not to cause us problems. The solution for this is a Gaussian function multiplied by a term that normalises it in our semi-infinite domain. This Gaussian will have a mean Vreset and a standard deviation σapprox „0which we will adjust according to our needs in future experiments. Thus, the equation (55) is: BPpVM, tq Bt“ ´BJpVM, tq BVM`JpVth, tqgpVMq(56) where gis the normalized Gaussian function that we have just defined. 53 Alejandro Barea Moreno ODE and PDE models used in neuroscience Lemma 5.3.1. The probability given by equation (56) is conserved when considering the boundary condition (54) at infinity. Proof: We can do the following: B BtżVth ´8 PpVM, tqdVM“żVth ´8 BPpVM, tq BtdVM,(57) due to the fact that the time is not the variable of integration and we assume that P is sufficiently regular. Then, substituting (56) into (57) we obtain: żVth ´8 BPpVM, tq BtdVM“ ´żVth ´8 BJpVM, tq BVM dVM`JpVth, tq “lim VMÑ´8 JpVM, tq`JpVth, tq´JpVth, tq “ lim VMÑ´8 JpVM, tq where we have used that: żVth ´8 gpVMqdVM“1. We know that the probability function being 0 at infinity imposes that the J function must also be 0, otherwise, seeing the definition of the flux of probability, the function P would not turn out to be 0. With this, we have the result. ˝ 5.3.2. θ-model Here things are easier. Let us remember that the θ-model was defined by equation: dθ dt “cp1´cospθqq` ab cp1`cospθqqpI´I1q. For now, consider a“b“c“1, I1“0. If we plug (48) into that equation, we obtain the following Langevin equation: dθt“f1pθt, µqdt `f2pθt, σqdBt(58) where: f1pθ, µq“p1´cospθqq`p1`cospθqqµ, f2pθ, σq“p1`cospθqqσ, which are both periodic functions. Therefore, the Fokker-Planck equation for the θmodel will 54 Alejandro Barea Moreno ODE and PDE models used in neuroscience have periodic coefficients and it is justifiable to add periodic boundary conditions: lim θÑπ´Ppθ, tq “ lim θÑ´π`Ppθ, tq,(59) lim θÑπ´Jpθ, tq “ lim θÑ´π`Jpθ, tq.(60) and, because of the periodicity of the coefficients and their derivatives, we have: Property 5.3.1. The second periodic boundary condition is equivalent to: lim θÑπ´BPpθ, tq Bθ“lim θÑ´π`BPpθ, tq Bθ. Then, we have the following: Lemma 5.3.2. The solution to the Fokker-Planck equation generated by (58) is conserved with the periodic boundary conditions (59) and (60). The proof is analogous to the one in Lemma 5.3.1. 5.4. Analysis of the solutions We analyze mathematically the solutions for the θmodel, as the case of the Leaky Integrate- and-fire includes a term JpVM, tqwhich is an evaluation of the solution we are trying to calculate and that ruins our efforts at some previous analysis. First of all, consider strong solutions to the problem: BPpθ, tq Bt“ ´Bpf1pθqPpθ, tqq Bθ`1 2B2pf2pθq2Ppθ, tqq Bθ2.(61) Lemma 5.4.1. The Fokker-Planck equation (61) with periodic boundary conditions (59) and (60) has an unique strong solution in the domain r´π, πsif the initial condition Ppθ, 0q “ P0pθq is analytic. Proof: Near t“0, the application of Cauchy-Kovalevska theorem (Theorem 0.6.5) tells us that the solution exists and is unique, because of the analyticity of the coefficients f1and f2, which are combinations of cosines. If then we consider that probability is conserved, the solution does not blow up and therefore we have global unique solutions. ˝ Calculating this solution can be very tedious as we would have to use analytic methods which would yield very complex solutions that are difficult to compute, so we would have to resort to the weak solutions. For this, we will use the following theorem found in [15]: 55 Alejandro Barea Moreno ODE and PDE models used in neuroscience Theorem 5.4.1 (Lions).Let VĂHbe two separable Hilbert spaces such that V is dense in H. Let M, α, T ą0, β ě0, F PL2p0, T ;Hq, u0PHand a:r0, T sˆVˆVÑVa family of bilinear applications, then if: |apt;u, vq|ďM∥u∥V∥v∥V,@u, v PV. apt;u, vq ě α∥v∥2 V´β∥v∥2 H,@vPV. tÑapt;u, vqis measurable @u, v PV. Then, the variational problem: $ ’ ’ & ’ ’ % Find uPL2p0, T;VqXC0pr0, T s;Hqsuch that d dt pu, vqH`apt;u, vq “ pFptq, vqH,@vPV, for almost all tP p0, T q up0q “ u0 has a unique solution. Now, we can discuss the weak solutions. Lemma 5.4.2. The Fokker-Planck equation (61) with periodic boundary conditions (59) and (60) has an unique weak solution in the domain r´π, πsif the initial condition Ppθ, 0q “ P0pθq is in L2pΩq. Proof: Let us start establishing the settings for the variational problem in which we will transform the PDE. Take H“L2pΩqwhere Ωis r´π, πs, and Vwill be the following subset of Sobolev space H1pΩq: V“␣vPH1pΩq|vp´πq “ vpπq( which is Hilbert because it is closed in H1pΩq. The norms in those two spaces are: ∥v∥H“żπ ´π |vpθq|2dθ. ∥v∥V“żπ ´π˜|vpθq|2`     dv dθ pθq     2¸dθ. Now, to obtain the bilinear form, we must choose a test function vPV, multiply and then integrate, to obtain (we will drop the dependencies for neatness): żBP Btv`żBpf1Pq Bθv´1 2żB2pf2Pq Bθ2v“0 56 Alejandro Barea Moreno ODE and PDE models used in neuroscience Fig. 35: Population activity over time of the θneuron model applied to a neuronal population. We see that the solution is periodic, and the waves will produce indefinitely, however, due to the numerical defects of the scheme we are using, the probability is lost little by little and therefore the wave is fading. The stationary solution is thus never attained in theory. This loss of probability is a function of the number of nodes used in time and in space. 63 Alejandro Barea Moreno ODE and PDE models used in neuroscience Conclusions Throughout this work we have seen a staggering variety of neuron models and some socalled neuronal population models for clusters of cells, but there is still a great deal to be seen, many mathematical artefacts of great calibre that, for the brevity of this work, must remain unexamined. Each chapter of this work acts as a node in the tree of computational neuroscience, because from each small Lemma, equation or concept, a large number of related concepts appear that have been or are being treated by some of the experts in the field, so I wanted to restrict myself to some of the most relevant aspects. This work does not pretend to be an exhaustive analysis of the subject and to exhaust it, since a field as extensive as this, which is growing day by day, is impossible to cover in so few pages. However, it is a brief introduction to some of the concepts on which current research of such calibre as the European Union’s flagship project, The Human Brain Project, is based, as well as such fundamental tools in the work of health researchers as The Virtual Brain. Fig. 36: Simulations of the brain done with The Virtual Brain [19]. We have also neglected concepts such as the detailed study of the waves that travel through the brain, the delayed integrodifferential equations that serve to explain communication over long distances or the models of memories based on neural networks, which are so fashionable today. A good synthesis of most of the knowledge in the field is found in [18] which is a paper that has influenced greatly the structure of this work. Still, I hope that the reader will be satisfied with my treatment of mathematics and neuroscience and that he or she will have been able to retain some of the ideas that were presented. 64 Alejandro Barea Moreno ODE and PDE models used in neuroscience Appendix Simulations code All of the code for the simulations have been developed by me using MATLAB. As I have used a lot of little scripts and including them here would prolong the work unnecessarily, I link you to my GitHub where I have uploaded everything: https :{{github.com{alejbarea{T F Gneuroscience.git Previous mathematical results We will need some results known to everyone that has taken part in courses on ODEs and stability theory. These will be stated in this section without proof. These results have been extracted from [17] and [16]. Definition 0.6.1 (Hyperbolic fixed point).Let F:RnÑRnbe a dynamical system with a fixed point p. Then if J, the jacobian matrix of Fat p, has no eigenvalue with zero real parts then pis hyperbolic. Theorem 0.6.1 (Local Existence and Uniqueness Theorem).Let ΩĂRN`1an open set and f:ΩÑRNso that fPC0pΩq X Liplocpy; Ωq. Under these conditions, for each pt0, y0q P Ω, there exists one and only one maximal solution to the Cauchy Problem: pCPq#y1“fpt, yq ypt0q “ y0 Lemma 0.6.1 (Gronwall’s Lemma).Let 0ăt0ăt1ă 8 and u, a PC0prt0, t1sq with aptq ě 0,@tP rt0, t1sand: uptq ď M`żt t0 apsqupsqds, @tP rt0, t1s. Then: uptq ď Meşt t0apsqds,@tP rt0, t1s. Theorem 0.6.2 (Sufficient condition for Lipschitzianity).Let ΩĂRN`1an open set and f“ pf1,¨¨¨ , fNq:ΩÑRN. a function whose partial derivatives Bfi Byjexists and are continuous @i, j “1,¨¨¨ , N then fPLiplocpy, Ωq. 65 Alejandro Barea Moreno ODE and PDE models used in neuroscience Theorem 0.6.3 (Stability by linearization).Let y1“fpyq be a differential system with fPC2pBpq, BpĂRN, N ě1a closed ball. Let y˚be a critical point of the system, i.e. fpy˚q “ 0. If the eigenvalues of the Jacobian matrix evaluated in y“y˚,pBfi Byjpy˚qqij have all negative real parts, then ϕ˚, the solution with an initial condition in y y0“y˚, is an exponentially asymptotically stable equilibrium of the system. If any of the eigenvalues has a positive real part, then the equilibrium is unstable. Theorem 0.6.4 (Poincaré-Bendixon).Let y1“fpyqbe an autonomous system with fP C1pDq, D ĂR2an open and connected non-empty set. Let y0PKĂDwith D a compact that γ`py0q Ă Kwith γ`py0qbeing the orbit of the system starting from t0. Then, if the set of limit points Λ`py0qdoes not have critical points of the system we have: Λ`py0qis a non-degenerated cyclic orbit. Either γpy0q “ Λ`py0qor γ`py0qapproaches Λ`py0qin spiral motion. In this case we call Γ`py0qalimit cycle. Theorem 0.6.5 (Cauchy-Kovalevska).A PDE of the form: BMu BxM 1“fpxi, u, Bm1`¨¨¨`mnu Bxm1 1¨¨¨Bxmn nq where m1`¨¨¨`mn“Mand boundary conditions for x1“x0 1: Bju Bxj 1“Ljpx2,¨¨¨ , xnq, j “0,¨¨¨ , M ´1. If f and Ljare analytic in a neighbourhood of upx0 1,¨¨¨ , x0 nqand Bm1`¨¨¨`mnu Bxm1 1¨¨¨Bxmn npx0 1,¨¨¨ , x0 nqthen the equation has an unique local solution in a neighbourhood of px0 1,¨¨¨ , x0 nq. Theorem 0.6.6 (Hölder).Let p, q P p1,8q such that: 1 p`1 q“1. And f, g :XÑ 8 are two measurable positive functions over X. Then: żfgdµ ď pżfpdµq1{ppżfqdµq1{q 66 Alejandro Barea Moreno ODE and PDE models used in neuroscience References [1] Levitan, I. B., &Kaczmarek, L. K. (2002). The Neuron: Cell and Molecular Biology. Oxford University Press. [2] Ermentrout, G. B., &Terman, D. H. (2010). Mathematical Foundations of Neuroscience. Springer Publishing. [3] Strogatz, S. H. (1998). Nonlinear Dynamics and Chaos. En International Workshop. https://doi.org/10.1142/3481. [4] Hodgkin, A. L., &Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of Physiology, 117(4), 500-544. https://doi.org/10.1113/jphysiol.1952.sp004764. [5] Izhikevich, E. M. (2006). Dynamical Systems in Neuroscience: The Geometry of Excitability And Bursting (1.a ed.). Mit Pr. [6] Hartman, P. (1960). A Lemma in the theory of structural stability of differential equations. Proceedings of the American Mathematical Society, 11(4), 610-620. https://doi.org/10.1090/s0002-9939-1960-0121542-7. [7] Troy, W. C. (1978). The bifurcation of periodic solutions in the Hodgkin- Huxley equations. Quarterly of Applied Mathematics, 36(1), 73-83. https://doi.org/10.1090/qam/472116. [8] Gerstner, W., Kistler, W. M., Naud, R., &Paninski, L. (2014). Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition. Cambridge University Press. [9] Novère, N. L. (2012). Computational Systems Neurobiology. Springer Publishing. [10] Øksendal, B. (2006). Stochastic Differential Equations: An Introduction with Applications. Springer Publishing. [11] Prados, A. (2022). Ampliación de Física Estadística. Universidad de Sevilla. [12] Rolls, T., &Deco, G. (2010). The Noisy Brain: Stochastic Dynamics as a Principle of Brain Function. OUP Oxford. [13] Herculano-Houzel, S. (2009). The human brain in numbers: a linearly scaled-up primate brain. Frontiers in Human Neuroscience, 3. https://doi.org/10.1113/jphysiol.1952.sp004764. [14] Gardiner, C. W., Gardiner, C. W., &G. (1985). Handbook of Stochastic Methods: For Physics, Chemistry and the Natural Sciences (2nd ed.). Springer. 67 Alejandro Barea Moreno ODE and PDE models used in neuroscience [15] Fernández Cara, E. (2023). Complementos de Modelización y Optimización numérica. Universidad de Sevilla. [16] Fernández Cara, E. (2023). Ecuaciones Diferenciales Ordinarias. Universidad de Sevilla. [17] Romero Enrique, J. M. (2021). Física Matemática. Universidad de Sevilla. [18] Deco, G., Jirsa, V. K., Robinson, P. N., Breakspear, M., &Friston, K. J. (2008). The Dynamic Brain: From Spiking Neurons to Neural Masses and Cortical Fields. PLOS Computational Biology, 4(8), e1000092. https://doi.org/10.1371/journal.pcbi.1000092 [19] The Virtual Brain. (s. f.). https://www.humanbrainproject.eu/en/medicine/the-virtual- brain/ 68