Computational simulation of a new system modelling ions electromigration through biological membranes

Background The interest in cell membrane has grown drastically for their important role as controllers of biological functions in health and illness. In fact most important physiological processes are intimately related to the transport ability of the membrane, such as cell adhesion, cell signaling and immune defense. Furthermore, ion migration is connected with life-threatening pathologies such as metastases and atherosclerosis. Consequently, a large amount of research is consecrated to this topic. To better understand cell membranes, more accurate models of ionic flux are required and also their computational simulations. Results This paper is presenting the numerical simulation of a more general system modelling ion migration through biological membranes. The model includes both the effects of biochemical reaction between ions and fixed charges. The model is a nonlinear coupled system. In the first we describe the mathematical model. To realize the numerical simulation of our model, we proceed by a finite element discretisation and then by choosing an appropriate resolution algorithm to the nonlinearities. Conclusions We give numerical simulations obtained for different popular models of enzymatic reaction which were compared to those obtained in literature on systems of ordinary differential equations. The results obtained show a complete agreement between the two modellings. Furthermore, various numerical experiments are presented to confirm the accuracy, efficiency and stability of the proposed method. In particular, we show that the scheme is unconditionally stable and second-order accurate in space.


Introduction
In this paper we consider a class of models of ions migration through biological membranes. Such migrations exist for most living cells and some biochemical processes. The motion of ions is supposed due to diffusion and to the effect of the electrical field. Furthermore ions can undergo reactions. So the ions concentrations satisfy the Nernst-Planck equations, including a kinetic reaction terms and the potential is given by Poisson equation. The model is given by: where Q T = ]0, T[ × , T = ]0, T[ × ∂ , T > 0; is an open regular set of R 2 which represent the biological cell and ∂ represent the cell membrane. For each i, C i is the concentration of the i species which has diffusion coefficient d i , mobility m i and valency z i . φ is the electrical potential, f is the fixed charges concentration http://www.tbiomed.com/content/10/1/51 and the F i are reaction terms. We suppose that F i depends continuously on the C j 's and φ, and that f is a bounded function. We suppose that d i is a positive constant for each i.
for all i = 1, . . . , Ns C i,0 ∈ L 2 ( ) and satisfy C i,0 ≥ 0 ( 2 ) Further information about the modelling of the problem and its mathematical analysis, can be found in [1].
In the biophysical literature, the early works on these models were interested in the stationary case of passive migration (i.e., without reaction); see [2,3]. In all these works two popular simplifications were considered, namely the Goldman hypothesis where the electrical field is supposed to be constant inside the membrane and the electroneutral hypothesis where the neutrality at each point of the membrane is assumed; see for example [3]. Mcgillivray [4] recognized that these models are the limit of the full equations when the ratio √ ε = λ l of the Debye length to the membrane thickness goes to, respectively, infinity or zero. Usually enzymes are held to biological membranes and ions undergo biochemical reactions when crossing the membrane. Valleton [5] did a general biophysical study of coupling of electromigration diffusion with biochemical reactions.
In this paper we present a numerical simulation of such systems, for a large class of reaction kinetics, including the usual biochemical kinetics as the Michaelis-Menton one (a mathematical analysis of the one dimensional and stationary case was done by [6] then we did the mathematical analysis of the multidimensional unsteady case [1]). This article is organized in the following way. The next section is devoted to finite element discretisation of the mathematical model. Then, we present applications, results and numerical experiments showing the accuracy, efficiency and stability of the proposed method. Finally, conclusions are drawn in the last section.

Finite element discretisation
In order to show the numerical formulation of the problem, let V = L 2 (0, T; H 1 ( )) be the space of approximate solutions and W = H 1 ( ) be the space of tests functions. Let W h be a finite element space of Lagrange P1 included in W and V h = L 2 (0, T; W h ) be the finite dimensional subspace of V. Now we introduce the function Moreover, we consider The Faedo-Galerkin formulation for the problem is given by, finding (3) http://www.tbiomed.com/content/10/1/51 Let (y j ) 1≤j≤m the mesh nodes, ( j ) 1≤j≤m the canonical basis of W h , we consider the following two sets of index Now let's consider an uniform subdivision of [0, T], we define a time step dt = T N , for N ≥ 1. We pose then: and ζ φ n the approximation of ζ φ (t n ), then we used an implicit scheme for the discretization of the time derivative. By the method of finite elements see [7], we arrive at the following formulation of the problem: Find the vectors nodal concentrations ξ Z n+1 Let T h the mesh generation of containing nel finite elements. For a finite element e k ∈ T h , let be T e k = {k 1 , k 2 , k 3 } where k 1 , k 2 , k 3 are the numbers of degrees of freedom of e k and N e k k l are interpolation functions. We We have a nonlinear term due to S i s ,ζ φ n+1 (ξ Z n+1 Ns ) , we have dealt with according to the model and thus to the expression of S i s ,ζ φ n+1 . http://www.tbiomed.com/content/10/1/51

Results and discussion
In this section we present three numerical applications of ions electro-migration through biological membranes. The models of the basic enzyme reaction, the suicide substrate reaction and the cooperative reaction, are numerically simulated.

Result 1 : Enzyme kinetics (basic enzyme reaction)
To understand where some of the more complicated reaction schemes come from, we consider a reaction that is catalyzed by an enzyme. Enzymes act as remarkably efficient catalysts (generally proteins), by accelerating the conversion of some other molecules called substrates into products via lowering the free energy of activation of the reaction, but they themselves remain unchanged by the reaction. Thus, they are important in the regulation of biological processes, for example as activators or inhibitors. One of the most basic enzymatic reactions, first suggested by Michaelis and Menten [8], implies a substrate S reacting with an enzyme E to form a complex SE which in turn is converted into a product P and the enzyme. This is represented schematically by here k 1 , k −1 and k 2 are constant parameters associated with the rates of reaction. We denote the concentrations of the reactants by Then the law of Mass Action applied to (4) leads to one equation for each reactant and hence a system of nonlinear equations. The usual approach to these equations is to assume that the initial stage of the complex C 4 , formation is very fast after which it is essentially at equilibrium, then we get C 4 in terms of C 2 ,

Algorithm of resolution
Before stating the resolution algorithm, we introduce the function Z h = Z i,h 1≤i≤3 defined by Moreover, we consider We used the following algorithm to calculate φ h and Z i,h then we calculate C i,h by using the reverse relation: • Loop over n At step n : where eps is the stopping criterion.

Numerical results
Here we present changes in substrate, product and enzyme concentrations. Shown in Figure 1 is the spatial distribution of substrate concentration through the cell at both initial and final time (t=0 and T=200 ms). Figure 2 presents the spatial distribution of the product concentration through the cell at both initial and final time (t=0 and T=200 ms).
In Figure 3, one sees the time evolution of the substrate and the product concentrations at the center of the cell. It can be seen that the substrate decrease curve is the mirror image of the product appearance curve. By observing the early times, it's obvious that the substrate loss and product appearance change speedily with time but as time goes on these rates diminish, to reach zero when all the substrate has been converted to product by the enzyme.
The Figure 4 shows, as predicted, the enzyme concentration remaining constant over time.
Two simplifications of these equations have been quite popular in the literature while computing membrane reactions, firstly the Goldman hypothesis where the electrical field is supposed to be constant inside the membrane and secondly considering a system of ordinary equations depending only on time and not the space. The added value of this work is not considering all of those simplifications which leads to a more realistic model and more accurate numerical results. Moreover, the results obtained are in agreement with the experimental results found in the literature [9].

Result 2 : Suicide substrate kinetics
An enzyme system of major experimental concern; see [10,11], is the mechanism-based inhibitor, or suicide substrate system, represented by Walsh et al. [12], where E, S and P stand for enzyme, substrate, and product, respectively; X and Y, enzymesubstrate intermediates; E i , inactivated enzyme; and the k 's are positive rate constants.
In this system, Y has a choice of one of two pathways, namely, to E + P with rate k 3 or to E i with rate k 4 . The ratio of these rates, k 3 /k 4 , is called the partition ratio and is denoted by r. Each of these pathways are supposed to be irreversible over the timescale of the reaction see [13]. S is known as a suicide substrate because it binds to the active site of an enzyme-like a substrate-but the enzyme converts it into an inhibitor which irreversibly inactivates the enzyme. Thereby, the enzyme 'commits suicide' . In this way, a suicide substrate can specifically target an enzyme for inactivation. Furthermore, suicide substrates are particularly useful in drug administration, as they are not noxious in their common form and only the designated enzyme can convert them to their inhibitor form. For example, suicide substrates have been subject of investigation for use in the treatment of depression (monoamine oxidase inhibitors, Seiler at al. [10]), epilepsy (brain GABA transaminase inhibitors, Walsh [11]), and some tumors (ornithine decarboxylase inhibitors, Seiler et al. [10]). Suicide substrate kinetics have been studied by Waley [13] and by Tatsunami et al. [14], who had interest in the factor which determined whether the substrate was exhausted before all the enzyme was inactivated. Waley proposed it was rμ, where μ is the ratio of the initial concentration of enzyme to that of substrate, namely, e 0 /s 0 . Tatsunami et al., on the other hand, found the determining factor to be (1 + r)μ. When (1 + r)μ > 1 the substrate is exhausted, while for (1 + r)μ < 1, all the enzyme is inactivated. When (1 + r)μ = 1, both occur. The interest is when e 0 /s 0 is not small, which was in effect assumed since both Waley [13] and Tatsunami et al. [14] used a quasi-steady state approximation. The validity decreases for increasing values of e 0 /s 0 . We denote the concentrations of the reactants by The law of mass action applied to (6) leads to one equation for each reactant and hence a system of nonlinear equations. We obtain the following system

Algorithm of resolution
Before stating the resolution algorithm, we introduce the function Z h = Z i,h 1≤i≤6 defined by . . , 6 http://www.tbiomed.com/content/10/1/51 Moreover, we consider We used the following algorithm to calculate φ h and Z i,h then we calculate C i,h by using the reverse relation: • Loop over n At step n :

Numerical results
Here we present changes in substrate, product, enzyme, inactivated enzyme and the intermediate concentrations (X and Y ). The time step of the simulation is dt = 10 −2 s. The data employed for the reaction parameters and initial concentrations were taken from Burke et al. [15]. Figure 5 plots the changes in the concentration distribution of substrate from initial time to final time. One can see that in final time the substrate was totally exhausted. This complete consumption of the substrate is in agreement with the prediction of Tatsunami et al. [14] as (1 + r)μ = 6 > 1.
In Figure 6, one sees the evolution of substrate concentration over time at the center of the specimen. Figure 7 shows the evolution of inactivated enzyme concentration over time at the center of the specimen. Figure 8 shows the numerical solutions for intermediate concentrations of X and Y.
In Figure 9, is represented the graphic of the evolution of substrate and product at the center of the cell, comparing that result with Figure 3 (Michaelis and Menten model), here the two plots are asymmetric which is logical as we know that an amount of the substrate instead of being converted to product, is forming the inactivated enzyme (inactivating the enzyme).
To illustrate, Figure 10 shows the decrease in the enzyme concentration unlike the Michaelis Menten model ( Figure 4); however, as the intermediate enzymes X and Y, vanish in few milliseconds, we see the loss in enzyme compensated by the production of the inactivated enzyme: the enzyme commits suicide.
To highlight the accuracy of these results, we compared them first with the numerical solutions and the approximate asymptotic solutions obtained both by Burke et al. [15], they considered a system of ordinary differential equations depending on time as they neglected the spatial aspect of the biochemical reaction, and they supposed that the electrical potential inside the membrane remains constant. For the numerical solutions Burke et al. [15] solved the system numerically, but for the approximate asymptotic solutions they non-dimensionalise the same system, and used asymptotic methods and a method detailed in Kevorkian and Cole [16]. Finally we compared our results with previous approximate methods of Tatsunami et al. [14] and Waley [17] which were based on a pseudo-steady state hypothesis. This comparison shows that the results described here are valid numerical solutions for the kinetics of suicide substrate system. The solution for the substrate and inactivated enzyme are more accurate than those of previous approximations [14,17], especially in small time, which is by definition ignored by any pseudo-steady-state approximate method. Furthermore, the method presented here is specially useful in estimating the intermediate (X and Y ) concentrations besides incorporating the spatial and the electro-migration aspects of the phenomena.

Result 3 : Cooperative phenomena
An enzyme reaction is said to be cooperative if a single enzyme molecule, after binding a substrate molecule at one site can then bind another substrate molecule at another site. Such phenomena are quite common in living organisms. Another interesting cooperative reaction is when an enzyme with several binding sites is such that the binding of one substrate molecule at one site can affect the activity of binding other substrate molecules at another site. This indirect interaction between distinct and specific binding sites is called allostery, and an enzyme displaying it, an allosteric enzyme. When a substrate that binds at one site increases the binding activity at another site then the substrate is called an activator, otherwise (if it decreases the activity) it's called an inhibitor.
As an example of cooperative phenomenon we consider the case when an enzyme has two binding sites. A model for this consists of an enzyme molecule E which binds a substrate molecule S to form a single bound substrate-enzyme complex X. This complex X not only breaks down to form a product P and the enzyme E again, it can also combine with another substrate molecule to form a dual bound substrate-enzyme complex Y. This Y complex breaks down to form the product P and the single bound complex X. The reaction mechanism is represented schematically by Here k s are rate constants. We denote the concentrations of the reactants by . http://www.tbiomed.com/content/10/1/51 Then the law of Mass Action applied to (8) leads to one equation for each reactant and hence the system of nonlinear reaction. We have then

Algorithm of resolution:
Before stating the resolution algorithm, we introduce the function Z h = Z i,h 1≤i≤5 defined by Moreover, we consider We used the following algorithm to calculate φ h and Z i,h then we calculate C i,h by using the reverse relation: • Loop over n At step n :

Numerical results
Here, we present the evolution of product and enzyme concentrations over time in both positive and negative cooperativity. We used the same constant parameters as the previous example for the diffusion coefficients of the species and the same initial conditions. The cell is represented by an ellipse with semi-major axis a=2 and semi-minor axis b=1.
To ensure cooperativity, the positive rate constants are chosen by the following reasoning: Suppose that the binding of the first substrate molecule is slow, but that with one site bound, binding of the second is fast (this is large cooperativity). This can be modeled by letting k 3 → ∞ and k 1 → 0 while keeping k 1 k 3 constant, in which case K 2 → 0 and K 1 → ∞ while K 2 K 1 is constant ( K 1 and K 2 were introduced as they appear in the http://www.tbiomed.com/content/10/1/51 expression of velocity reaction, this is discussed in details in the book by Keener and Sneyd [18]), where An enzyme can also exhibit negative cooperativity, in which the binding of the first substrate molecule decreases the rate of subsequent binding. This can be modeled by decreasing k 3 . We used for positive cooperativity K 1 = 1000, K 2 = 0.001 and K 1 = 0.5, K 2 = 100 for negative cooperativity (this values were taken from [18]).    In Figure 11, one sees that in positive cooperative reaction, the product concentration is characterized by an "S-shaped" sigmoidal curve, which is different from other enzyme reaction that exhibits a curve that tends to be hyperbolic. This results from cooperative effects; in which the enzyme can bind more than one substrate molecule, but the binding of one substrate molecule affects the binding of subsequent one.
In Figure 12, we plot the evolution of the enzyme concentration for both extreme positive cooperativity and negative cooperativity.

The convergence test
The rate of convergence of the scheme is difficult to prove analytically. However, numerical experimentation suggests that the scheme is second-order accurate in space. A quantitative estimate of the convergence error was obtained by performing a number of simulations for the same initial condition on a set of increasingly finer space meshes and time steps. The initial conditions are constants. Let T h the mesh generation of , and h(T h ) = max{diam(e k )|e k ∈ T h }, we take h = 0.3, h = 0.1 and h = 0.05. For each mesh we integrate to time T with dt = hT 16 . Note that as we refine the space step we also refine the time step. The error of the numerical solution was defined as

The convergence of the basic enzyme reaction
In Table 1 is presented the error of convergence for different mesh sizes in the case of basic enzyme reaction.  in subfigures (a), (b) and (c). The time steps are given below each subfigure.

The convergence of the suicide substrate reaction
In Table 2 is presented the error of convergence for different mesh sizes in the case of suicide substrate reaction.

The convergence of the cooperative reaction
In Table 3 is presented the error of convergence for different mesh sizes in the case of cooperative reaction.

Stability and accuracy tests
Now, let us give some information about the numerical stability of our algorithms. We perform a numerical experiment with different time step dt, dt 2 and dt 4 . These results suggest that the scheme is indeed unconditionally stable as the solutions are quasi the same for different time steps. To illustrate, we chose to represent the product concentration.

Stability of the basic enzyme reaction
In Figure 13, We display snapshots of the product concentration at time T=0.5 with three different time steps dt=0.00025, dt=0.0005 and dt=0.001. We can see that the results are quasi the same at the final time T.

Stability of the suicide substrate reaction
In Figure 14, we display snapshots of the product concentration at time T=4 with three different time steps dt=0.0025, dt=0.005 and dt=0.01. We can see that the results are quasi the same at the final time T. Figure 15 shows snapshots of the product concentration at time T=4 with three different time steps dt=0.0025, dt=0.005 and dt=0.01. We can see that the results are quasi the same at the final time T.  in subfigures (a), (b) and (c). The time steps are given below each subfigure. http://www.tbiomed.com/content/10/1/51

Conclusion
In this paper, a new model simulating ions electro-migration through biological membranes is proposed by using a more general mathematical model and a numerical technique based on the finite element method. The results presented here demonstrate that the model's behavior agrees with the behavior of biochemical reactions as it's consistent with the physical interpretation of the phenomena. Moreover, after comparison we can observe a complete consistency with literature findings [9,[13][14][15]18]. A variety of numerical experiments were presented to confirm the accuracy, efficiency, and stability of the proposed method. In particular, the scheme was shown to be unconditionally stable and second-order accurate in space.