Numerical algorithm for simulations of Brownian particles in an ionic channel with fixed-density boundary conditions

Numerical algorithm for simulations of Brownian particles in an ionic channel with fixed-density boundary conditions Oscar Escolano, Eric Planas (Dated: May 31, 2019) We learn a new programming language, Fortran, in order to use it to simulate a system of many particles that follow Brownian dynamics. In particular, we have modelled an ionic channel under the action of different potentials and we have analysed the ion density profile during both transient and steady states and compared it with the theoretical solution of the Fokker-Planck equation. The algorithm includes the temporal evolution of each individual particle through the integration of an stochastic differential equation and the implementation of fixed-density boundary conditions. I. INTRODUCTION This project has two main objectives. The first one is to learn Fortran, a new programming language for us, and the second one is to apply it to do our first simulation of a complex system. Specifically, we have simulated an ionic channel. This kind of structures span the cell membrane and are very relevant in cellullar physiology, since thanks to the fast ion flows through them, they allow the neurons to depolarize and transmit the action potential. It is interesting, then, to develop and algorithm that allows us to study this system under different membrane potentials and ion concentrations in and out of the cell in order to better understand its dynamics. II. THEORETICAL FRAMEWORK Our approach to simulate the system will be based in the proposed one by our tutor Laureano Ramirez de la Piscina in his work [1]. Hence, we will consider a set of independent particles that move along the channel. It will be done, however, a translation of the problem to one-dimension provided that the system has a paraxial symmetry and the electric field that appear is parallel to the channel axis. We will also consider that the dynamics of each particle in a channel of length Lwill follow the Langevin equation. Considering a small Reynolds number, it can be assumed the overdamped limit and use the first order equation γdx dt =−dV (x) dx +ξ(t) (1) Where ξ(t) is a Gaussian noise with zero mean and a standard deviation √2kBTγ that will simulate the thermal fluctuations and other stochastic forces that may act on the particles, γis the friction coefficient and Vis a time-independent potential that will simulate the membrane potential. Regarding the boundary conditions, we will consider that at both sides of the channel there are two particle reservoirs with a given fixed values of densities, since the flow of ions in and out of the channel will suppose minimum changes in the values of the concentration of this ion in and out of the cell. Once we obtain the results of the simulations applying the above-mentioned model, we will compare them with the Fokker-Planck equation for concentrations as a theoretical reference ∂ρ(x, t) ∂t =−∂ ∂xf(x)ρ(x, t) + D∂2 ∂x2ρ(x, t) (2) where f(x) and Dare given by f(x) = 1 γ dV (x) dx , D =kBT γ(3) III. NUMERICAL ALGORITHM In order to simulate the system, we have written a code that will take into account different items such as the temporal evolution of the particles, the implementation of the boundary conditions and the density calculations that will be explained below. All of them can be checked in the code in Appendix A. A. Time-discretized Langevin dynamics In order to implement the simulation of the particle motion in an ionic channel, many aspects have to be taken into account in parallel. Firstly, it must be considered the time evolution of the system. As stated above, we will use an Euler algorithm to integrate the overdamped Langevin equation. For this, a temporal step ∆thas been set, and the particle positions at the end of each time-step will be calculated using the previous position, obeying x(t+ ∆t) = x(t) + f[x(t)]∆t+1 γ√∆t ξ(t) (4) Where ξis a Gaussian random number with zero mean and standard deviation √2kBTγ. B. Boundary conditions As mentioned, fixed-density boundary conditions will be considered in order to study the dynamics at the channel boundaries. We will consider that at the end of each time step, the particles that have moved to a position out of the limits, that is if xis smaller than 0 or larger than L, will be removed from the system. That means that the boundaries will have no memory about the previous time iterations. To implement that, every time that a particle exits, its information must be replaced by the data of the last particle in the data vector that is still located in the channel. On the other hand, at any time step it must be also considered the possible entrance of particles in the channels. For the particles entrance, not only the rate at which a particle can enter has to be studied but also the position where it may enter. The mean number of particles that enter the channel in a time ∆twill be hNi=ρ√D∆t q(a) (5) Where a and the function q(x) are given by a=−fr∆t 4D(6) q(x) = −x erfc(x) + 1 √πe−x2(7) Note that the rate of entering particles will depend on the potential value at the boundary and the densities of the reservoirs. The number of particles that enter the channel each time should be computed from a Poisson distribution with the calculated mean number hNisince each particle arrival is an independent event with respect to the arrivals of other particles. However, with a small time step and not very large densities it can be seen that hNi  1. In this case the probability of appearing one single particle in a time ∆tcan be written as ρ(1) ≈ hNi. Then, in the algorithm it will only be necessary to generate a uniform random number χ∈(0,1) for each boundary and compare it with ρ(1). If χ≤ρ(1) then a particle will enter. Otherwise, no particle enters. In the case that in a certain time-step a particle enters, it should be decided at which position it is spawned. They cannot be placed just at the edge of the channels since due to the Brownian motion in very few iterations they would exit the channel again. To do this, the cumulative distribution function for the positions of entering particles has been used: F(x)=1−q(x−f∆t √4D∆t) q(a)(8) Then generating a uniform random number for each new particle, we get the value of Fand we obtain the desired position, which will be given by x=F−1(x) (9) In our case, at the beginning of the program, using Newton method, a library with many values of x and F(x) is created and whenever is needed the inverse function is computed via interpolation. There is, however, one exception, for large x, where to compute the corresponding x the following approximation is used the following expression: x≈f∆t+q−4γkBT∆tln (√π q(a)(1 −χ)) (10) Note that taking into account the possible drifts at the boundaries, we expect that the rate of entering and leaving particles in the boundaries at steady-state is constant, so the ion density near the boundaries should be the same than the densities in the corresponding reservoirs. C. Density calculation To calculate the particle density profile in the channel at a given time, we have made an algorithm that divides the channel in 500 bins and the number of particles in each bin at each time step is counted. As stated, our code uses several approximations provided that it is being considered that the densities both in the boundaries and inside the channel are not large and the time step is small enough. With these parameters, however, the measured densities are expected to present large fluctuations. In order to reduce the dispersion in the results it is convenient to average the results by running many independent realizations of the system reducing dispersion by 1/√N, being N the number of independent realizations. This can be done without problems since the particles are independent and noninteracting. D. Random Gaussian distribution generation To generate the Gaussian random variables we use the Box-M¨uller transform on a uniform random variable. To generate the uniform random variable we use a linear congruential generator for a seed of 8 bytes of memory. At each realization the seed is updated, giving new pseudo-random numbers. [3] 2 0 0.5 1 1.5 2 2.5 3 3.5 4 Distance(nm) 5 10 15 Ion density(particles/nm) Simulated Theoretical FIG. 1. Particle concentration vs position in the steady state for a null external potential and ρ1=ρ2= 10. Black line: numerical solution of Fokker-Planck equation under these conditions. Red line: Average over 10000 independent realizations in a temporal window T= 10 µs. It can be seen that the system almost reaches steadystate. IV. RESULTS In order to study the temporal evolution and the steady-state of this system, we have considered several cases with different potentials and densities at the boundaries. These different conditions will be set by changing the corresponding parameters in the code. In every case the same initial conditions have been set. That is, at the beginning of the simulations, there will be 50 uniformly distributed particles inside the channel. That means the initial density profile will be constant with ρ= 12,5. We have performed several simulations of Brownian particles in a channel of the following characteristics: The length has been set to L= 4, γ= 1000 , and the chosen temperature has been kBT= 25. We have set the system of units based on length 1 nm, energy 1 meV and time 0,1µs. We have also run the simulation during enough time for the system to be close to its steady-state. As stated before, all the simulation results will be compared with the numerical solution of Fokker- Planck equation, which will be painted in black. In Appendix B it is shown the code for its numerical solving. We first consider that there is no membrane potential in order to study the particle motion only taking into account the Brownian dynamics. In the first run of the simulation, we have considered that there is no drift (V= 0) and the ion concentrations at both 0 0.5 1 1.5 2 2.5 3 3.5 4 Distance(nm) 0 2 4 6 8 10 12 14 Ion density(particles/nm) Simulated(t=3us) Simulated(t=6us) Simulated(t=9us) Simulated(t=12us) Simulated(t=15us) Theoretical FIG. 2. Temporal evolution of particle concentration vs position until reaching steady-state for a null external potential and ρ1= 10 and ρ2= 1. Black line: numerical solution of Fokker-Planck equation under these conditions. Coloured lines: Average over 10000 independent realizations at different times. sides of the channel is the same (ρ1=ρ2= 10) (Fig. 1.). It can be appreciated that the result close to the steady state is a constant density profile that converges to the theoretical solution. It can be appreciated, however, that under these conditions the noise in the density calculation is quite strong and more realizations may be needed to reduce it. We have also run a simulation maintaining the null drift but changing the particle concentration at the boundaries ρ1= 10 and ρ2= 1. The results are shown in Fig. 2. where it can be seen that the particle density at the steady-state describes a linear profile, which also converge to the theoretical solution. In the following simulations there have been also considered several different boundary conditions but now a constant drift has been added (V= (qφ/L)x), simulating an standard membrane potential. In circumstances where the particle density at the boundaries is the same, a constant flow of particles is expected, so the density profile would be the same as in the case for null potential. We have decided, however, to simulate a more realistic situation where the concentrations in and out of the cell are different. This would correspond to a situation where the channel is open and the ion flow through it is governed by the electrochemical gradient. In Fig. 3 it is shown the temporal evolution of the system when ρ1= 10 and ρ2= 1. In Fig 4 these conditions have been exchanged (ρ1= 1 and ρ2= 10). In both cases it can be observed that at every time the solutions of the simulation converge to the theoretical solution. 3 0 0.5 1 1.5 2 2.5 3 3.5 4 Distance(nm) 0 5 10 15 Ion density(particles/nm) Simulated(t=2us) Simulated(t=4us) Simulated(t=6us) Simulated(t=8us) Simulated(t=10us) Theoretical FIG. 3. Temporal evolution of particle concentration vs position until reaching steady-state for a linear external potential with qφ = 8kBTand ρ1= 10 and ρ2= 1. Black line: numerical solution of Fokker-Planck equation under these conditions. Coloured lines: Average over 10000 independent realizations at different times. 0 0.5 1 1.5 2 2.5 3 3.5 4 Distance(nm) 0 5 10 15 Ion density(particles/nm) Simulated(t=2us) Simulated(t=4us) Simulated(t=10us) Theoretical FIG. 4. Temporal evolution of particle concentration vs position until reaching steady-state for a linear external potential with qφ = 8kBTand ρ1= 1 and ρ2= 10. Black line: numerical solution of Fokker-Planck equation under these conditions. Coloured lines: Average over 10000 independent realizations at different times. Finally, we have considered how to model the shape of a gate in the ionic channel. This problem is intrinsically a third-dimensional problem, since it involves the three-dimensional structure of the proteins that form the channel. Nevertheless, an equivalent projection to a one-dimensional problem can be done by adding a barrier potential. The shape of the barrier can be very diverse, but we, following 0 0.5 1 1.5 2 2.5 3 3.5 4 Distance(nm) 0 10 20 30 40 50 60 70 80 90 Ion density(particles/nm) Simulated(t=3us) Simulated(t=6us) Simulated(t=9us) Simulated(t=12us) Simulated(t=15us) Theoretical FIG. 5. Particle concentration vs position during the transient with a linear external potential plus a gaussian barrier specified in Eq.(11) and boundary densities ρ1=ρ2= 10. Black line: numerical solution of Fokker- Planck equation under these conditions. Coloured lines: Average over 10000 independent realizations at different times. It can be seen how particles accumulate in the barrier and can barely pass, as if the channel gate was closed. Ref [2], have opted for adding a Gaussian barrier to a constant drift. That is V(x) = qφ Lx+Vbe− (x−xb)2 2l2(11) Where Vbis the height, lis the width and xbis the centre position of the barrier. We have taken the values Vb= 8kBT,l=L/16 and xb=L/2. In Fig. 5 the temporal evolution of this system is shown. It is observed again how the simulations agree with the Fokker-Planck simulations. V. CONCLUSIONS We consider that the main goals of these project have been achieved: We are now fluent in programming in Fortran and we have managed to obtain the expected results in our simulations. We have been able to translate the problem into a numerical algorithm that, as checked, works well and fulfills the theoretical solutions of the Fokker-Planck equations. In order to go further in this project, it would also have been interesting to program the full three-dimensional problem. Nevertheless, this implies the implementation of more complex boundary conditions in yand zdirections, for which we would need more time. However, as the obtained results in the one-dimensional problem have been satisfying, we consider that having reached that point the presented project is complete by itself. 4 References [1] L Ram´ırez-Piscina. “Fixed-density boundary conditions in overdamped Langevin simulations of diffusion in channels”. In: Physical Review E 98 (Apr. 2018). doi:10.1103/PhysRevE.98.013302. [2] L Ram´ırez-Piscina and J Sancho. “Physical properties of voltage gated pores”. In: The European Physical Journal B 91 (Apr. 2018), p. 10. doi:10.1140/ epjb/e2017-80569-5. [3] W.T.Vetterling W.H.Press S.A.Teulosky and B.P.Flannery. Numerical Recipes: The Art of Scientific Computing. New York: Cambridge University Press, (2007). 5 Appendix A Simulation code program simulation ! Variable declaration : integer*8 :: i, nm , ns0 , nsF , ne0 , neF , nrealitzacions , ne0n , neLn integer*8 :: gg , itt , it1 , it2 , it3 , it4 , ipos real*8 :: N1 , N2 , rho1 , rho2 , D, kbt , gamma , dt , a1 , a2 , f1 , f2 real*8 :: Length , pot , derpot , qfun , arg , qa1 , qa2 , Fx , pi , Fxv real*8 :: suma2 , deltax real*8 :: x0 , y, xn , xre , newton , m, xentrada , dum , a, b, c, ns0n real*8 , dimension (0:9999) :: vFxL , vFy0 , vFyL , vFx0 integer*4 :: seed parameter (numbin=1000) parameter ( nparticmax =1000) parameter ( Length =4.0) parameter (dt =5.0E -4) real*8 , dimension ( nparticmax ) :: xpos real*8 , dimension (0: numbin -1) :: density , density1 , density2 real*8 , dimension (0: numbin -1) :: density3 , density4 , density5 real*8 , dimension (0: numbin -1) :: densityt1 , densityt2 , densityt3 real*8 , dimension (0: numbin -1) :: densityt4 , densityt5 dimension nbin (0: numbin ) real*8 , dimension ( nparticmax ) :: soroll real*8 :: U1 , U2 , r8_uniform_01 , r8_normal_ab , nsLn parameter ( kbt =25. d0 , gamma =1000. d0 , D= kbt / gamma ) parameter (rho1=10.d0, rho2=10.d0) parameter (itt =300000 , it1= itt/5, it2 =2* itt /5) parameter (it3=3*itt/5, it4=4*itt/5) common pi , seed pi =4. d0* atan(1.0 _8) seed =933543 print *,’ Initial ␣ number ␣ of␣ particles ␣=␣’, nparticulas deltax=Length/float(numbin)!Bin length ! Boundary conditions f1 =- derpot (0. d0 )/ gamma f2 =- derpot ( Length )/ gamma print *,’pot ␣0␣=␣ ’, pot (0. d0) print *,’pot ␣2␣=␣ ’, pot (2. d0) print *,’pot ␣L␣=␣ ’, pot (4. d0) print *, ’derpot0 ␣=␣’, derpot (0. d0 ) print *, ’derpotL ␣=␣’, derpot ( Length ) a1=-f1* sqrt( dt /(4. d0*D)) N1= rho1 *sqrt(D* dt )* qfun ( a1 ) a2= f2* sqrt(dt /(4. d0*D)) N2= rho2 *sqrt(D* dt )* qfun ( a2 ) print *, ’N1 ’,N1 print *, ’N2 ’,N2 1 ! Creation of a LIBRARY with the yequispaced nodes of the distribution ! function ( vFy0 ) and the corresponding x nodes ( xre ) at position 0: qa1 =- a1 * erfc ( a1 )+ exp(-( a1 **2))/ sqrt(pi) !q(a) at 0 x0 =5.3E -5 ! Initial guess y=0. d0 !We know that the first components of vFx0 and vFy0 are (0 ,0) !to guarantee the convergence of Newton ’s method vFx0 (0)=0. d0 vFy0 (0)=0. d0 do i =1 ,9999 y=y+1. d0 /10000. d0 !y equispaced xn = newton (y ,x0 , qa1 ) !It is the argument of qfun ( see eq .7) xre = sqrt(4. d0 *D* dt )* xn +f1 *dt ! Calculation of the real x ! Storage of the found values at this iteration vFx0 (i )= xre vFy0 (i)=y x0= xn ! The initial guess of the following iteration will be the previous ! found value with Newton ’s method . enddo ! Writing of the found vectors in a .dat in 2 rows open (32 , file=’ partentr0n .dat ’,status=’unknown ’) do i =0 ,9999 write (32 ,*) vFx0(i),vFy0 (i) enddo close (32) ! Creation of a LIBRARY with the equispaced nodes of the ! distribution function ( vFy0 ) and the corresponding x nodes !( xre ) at position Length . qa2 =- a2 * erfc ( a2 )+ exp(-( a2 **2))/ sqrt(pi) !q(a) at Length x0 =5.3E -5 ! Initial initial guess y=0. d0 !We know that the first components of vFxL and vFyL are (4 ,0) to ! guarantee the convergence of Newton ’s method vFxL (0)=4. d0 vFyL (0)=0. d0 do i =1 ,9999 y=y+1. d0 /10000. d0 xn = newton (y ,x0 , qa2 ) xre = sqrt(4. d0 *D* dt )* xn -f2 *dt ! Changes sign because we are evaluating at !Length vFxL (i )= Length - xre vFyL (i)=y x0= xn enddo ! Writing of the found vectors in a .dat in 2 rows open (36 , file=’ partentrFn .dat ’,status=’unknown ’) write (36 ,*) vFxL write (36,*) write (36 ,*) vFyL close (36) 2 do nrealitzacions=1,5000 ! Number of realizations of the simulation !Counters ne0n =0 neLn =0 ns0n =0 nsLn =0 ! Initial particle positions nparticulas =50 ! Places the particles in the indicated distribution !We use an equispaced initial distribution call coloca ( nparticulas , Length , xpos ) ! BEGINNING OF THE SIMULATION do i =1 , itt ! Number of temporal iterations ! Generation of the Gaussian noise using Box - Muller algorithm if (nparticulas >0) then do j=1, nparticulas soroll (j )= r8_normal_ab (0. d0 , sqrt(2. d0 * kbt * gamma ), seed ) enddo else print *,’Uep ␣no␣hi␣ha␣ particules !␣i=’, i endif ! Integration of the stochastic differential equation using Euler ’s !method if (nparticulas >0) then do j=1, nparticulas xpos (j)= xpos(j)+ sqrt( dt )*( soroll (j)/ gamma )- & dt * derpot ( xpos (j ))/ gamma enddo endif ! Generation of two random numbers in a uniform distribution [0 ,1] U1 = r8_uniform_01 (seed ) U2 = r8_uniform_01 (seed ) ! Particle ENTRY from x =0 if (U1 <N1) then nparticulas=nparticulas+1 ne0 =ne0 +1 ne0n=ne0n+1 !Each time that enters a particle , this counter increases by 1 3 U1 = r8_uniform_01 (seed ) !Xi. Different number than before to avoid !correlation ! Calculation of the position where it enters : if (U1 <=0.9999) then ! This is the maximum possible value of the !library nm = floor ( U1 *10000) ! 10000 because the Newton library has 10000 ! components m=( vFx0 (nm +1) - vFx0 (nm ))/( vFy0 ( nm +1) - vFy0 ( nm )) ! Slope ^ -1 xentrada = vFx0 ( nm )+ m*( U1 - vFy0 (nm )) ! Position where the particle !enters xpos ( nparticulas )= xentrada else ! Approximation dum =log(2. d0 *sqrt(pi )* qfun (a1 )*(1. d0 -U1 )) xentrada =f1 *dt+ sqrt( -4. d0 * gamma * kbt * dt *dum ) xpos ( nparticulas )= xentrada endif endif ! Particle ENTRY from x =4: if (U2 <N2) then nparticulas=nparticulas+1 neF =neF +1 ! Each time that enters a particle , this counter increases ! by 1 neLn=neLn+1 U2 = r8_uniform_01 (seed ) !Xi. Different number than before to avoid !correlation if (U2 <=0.9999) then ! This is the maximum possible value of the !library nm = floor (U2 *10000. d0)! 10000 because the Newton library has 10000 ! components m=( vFxL (nm +1) - vFxL (nm ))/( vFyL ( nm +1) - vFyL ( nm )) xentrada = vFxL ( nm )+ m*( U2 - vFyL (nm )) xpos ( nparticulas )= xentrada else dum =log(2. d0 *sqrt(pi )* qfun (a2 )*(1. d0 -U2 )) xentrada =Length -( -f2 *dt+ sqrt( -4. d0 * gamma * kbt * dt *dum )) xpos ( nparticulas )= xentrada endif endif ! Particle EXIT through the boundaries : if (nparticulas >0) then do j=1, nparticulas 4 endif enddo close (78) end program function pot (x) ! Different potentials used : real*8 :: x, pot , kbt , Length kbt =25. d0 Length =4. d0 pot =8. d0 *kbt / Length *x !pot = -50. d0*x +200. d0* exp (-(x -2. d0 )**2/(2. d0 *(0.25**2))) !pot =0. d0 return end ! Calculation of the numerical derivative of the potential : function derpot(x) real*8 :: x, derpot , dx , pot !dx =1E -8 derpot =( pot(x+dx)-pot (x-dx ))/(2* dx) ! Notice that it can fail when we ! have a non - derivable potential profile return end 11