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

Abstract

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.

Background

Cell membrane is the biological membrane separating the intracellular environment from the extracellular one. The cell membrane surrounds all cells and it’s selectively permeable, permitting the free passage of some substances and restricting the passage of others, thus controlling the flux of substances in and out of the cell. All diseases are problems of regulating the passage of materials at the level of the cell. Consequently, to understand the cause of a disease, we need to understand the alterations that take place at the cellular level. Thanks to mathematical modeling, cell biology phenomena may be expressed by ordinary differential equations or systems of partial differential equations.

An important class of models in cell biology, is the class modeling ion transport through biological membranes. This transport phenomena occurs in most living cells and some biochemical processes. The first models in the literature included one ordinary differential equation for each ion concentration. All this models were based on the implicit assumptions that chemical concentrations are uniform in space. This assumption is reasonable when the region of space where the reaction occurs is confined and very small. Also this models assumed that the electric field is constant inside the membrane.

Despite that the constant electrical field assumption has the advantage of leading to a simple mathematical analysis, all cells maintain a difference in electrical charge across their membrane. This difference in charge give rise to a voltage difference, or electrical potential. Furthermore, there are numerous situations in which chemical concentrations are nonuniform in space. In this sense, we need to establish and compute more accurate mathematical models of ions electromigration through biological membranes.

In this paper, we present the numerical simulation of a more realistic model of ions electromigration through biological membranes. This model is more general than those in literature of membrane transport as it extends them in four topics: 1) it’s a multidimensional model, 2) it doesn’t rely on the constant electrical field assumption, 3) it consider both the temporal and spatial dependence, 4) it includes different reaction kinetics terms.

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:

(1)

where Q _T =]0,T[×Ω, , T>0; Ω is an open regular set of 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 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.

(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 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²(0,T;H¹(Ω)) be the space of approximate solutions and W = H¹(Ω) be the space of tests functions. Let W^h be a finite element space of Lagrange P 1 included in W and V^h=L²(0,T;W^h) be the finite dimensional subspace of V. Now we introduce the function defined by

Moreover, we consider

The Faedo-Galerkin formulation for the problem is given by, finding Z_i,h∈V _h for i = 1,…,N _s and ϕ _h ∈V _h such that ϕ _h =0 in ∂ Ω:

(3)

where and are the projections of C_i,0, ϕ₀ on W _h .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

We represent the solutions as for i _s =1,…,N _s and with α _j (t)=0 for j∈Y₀. Then we set for i _s =1,…,N _s and ζ _ϕ =ζ _ϕ (t)=(α _j (t))_1≤j≤m.

Now let’s consider an uniform subdivision of [0,T], we define a time step for N≥1. We pose then:

Let’s note the approximation of and 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 for every 1≤i _s ≤N _s and nodal potential for j∈Y₀ such that

Let T _h the mesh generation of Ω containing nel finite elements. For a finite element e _k ∈T _h , let be where k₁, k₂, k₃ are the numbers of degrees of freedom of e _k and are interpolation functions. We have where s,l∈{1, 2, 3} if and only if Where , , for i _s =1,…,N _s , Q = (q _ij )_1≤i,j≤m, , β = (β _i )_1≤i≤m(X_Y×Y means the extracted matrix from X by keeping lines and columns with numbers belonging to Y)

where

where

where

where

where

where

Finally, finding for every 1≤i _s ≤N _s and for j∈Y₀ such that

We have a nonlinear term due to , we have dealt with according to the model and thus to the expression of .

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

(4)

here k₁, k₋₁ and k₂ 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₄, formation is very fast after which it is essentially at equilibrium, then we get C₄ in terms of C₂,

The basic enzyme reaction model becomes then

(5)

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:

● Initialize: for i = 1,2,3,

● Loop over n

At step n:

● Calculate solution of

● Calculate by:

● Calculate solution of

● Calculate solution of:

Loop over k untill

where eps is the stopping criterion.

● Calculate solution of

Numerical results

Here we present changes in substrate, product and enzyme concentrations. The cell is represented by an ellipse with semi-major axis a=2 and semi-minor axis b=1. The diffusion coefficient of the ions are d₁=10⁻³m².s⁻¹, d₂=2.10⁻³m².s⁻¹ and d₃=5.10⁻³m².s⁻¹. The constants of reaction are . The charge number of the ions are z₁=1, z₂=0 and z₃=1. The electric charge density is f = 0.1C. The initial conditions are C_1,0=1μ M, C_2,0=800 μ M, C_3,0=0 and ϕ₀=−80 m V; the stopping criterion is e p s = 10⁻⁴. The time step of the simulation is d t = 10⁻³s.

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 1

The spatial distribution of the substrate concentration through the cell at initial time t=0 and final time T=0.2 s.

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).

Figure 2

The spatial distribution of the product concentration through the cell at initial time t=0 and final time T=0.2 s.

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.

Figure 3

Evolution in time of substrate and product concentrations at the center of the cell.

The Figure 4 shows, as predicted, the enzyme concentration remaining constant over time.

Figure 4

Temporal evolution of the enzyme concentration at the center of the cell.

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],

(6)

where E, S and P stand for enzyme, substrate, and product, respectively; X and Y, enzyme-substrate 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₃ or to E _i with rate k₄. The ratio of these rates, k₃/k₄, 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₀/s₀. 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₀/s₀ 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₀/s₀. 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

(7)

Algorithm of resolution

Before stating the resolution algorithm, we introduce the function Z _h =(Z_i,h)_1≤i≤6 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:

● Initialize for i = 1,…,6

● Loop over n

At step n:

● Calculate solution of

● Calculate by:

● Calculate solutions of:initialize for i = 1,…,6,Loop overk untill

Numerical results

Here we present changes in substrate, product, enzyme, inactivated enzyme and the intermediate concentrations (X and Y). The cell is represented by an ellipse with semi-major axis a=2 and semi-minor axis b=1. The diffusion coefficients of the ions are d₁=10⁻³m².s⁻¹, d₂=2.10⁻³m².s⁻¹, d₃=5.10⁻³m².s⁻¹, d₄=10⁻³m². s⁻¹, d₅=2.10⁻³m².s⁻¹, d₆=4.10⁻⁶m².s⁻¹, the reaction parameters are k₁=2 s⁻¹, k₋₁=4 s⁻¹, k₂=12 s⁻¹, k₃=10 s⁻¹ and k₄=2 s⁻¹. The charge number of the ions are z₁=1, z₂=0, z₃=1, z₄=1, z₅=1 and z₆=0. The electric charge density is f = 0.1C. The initial concentrations are e₀=0.5μ M and s₀=0.5μ M;and ϕ₀=−80 m V. The time step of the simulation is d t = 10⁻²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.

Figure 5

Substrate concentration from initial time to final time.

In Figure 6, one sees the evolution of substrate concentration over time at the center of the specimen.

Figure 6

Evolution of substrate concentration over time at the center of the cell.

Figure 7 shows the evolution of inactivated enzyme concentration over time at the center of the specimen.

Figure 7

Evolution of inactivated enzyme concentration over time at the center of the cell.

Figure 8 shows the numerical solutions for intermediate concentrations of X and Y.

Figure 8

Evolution of the intermediate concentrations X and Y over time at the center of the cell.

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).

Figure 9

Evolution of substrate and product concentrations over time.

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.

Figure 10

Evolution of enzyme and inactivated enzyme over time at the center of the cell.

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

(8)

Here k^′s are rate constants. We denote the concentrations of the reactants by

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

(9)

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:

● Initialize for i = 1,…,5● 

● Loop over n

● At step n:

● Calculate solution of

● Calculate by:

● Calculate solutions of:initialize for i = 1,…,5,

Loop over k untill

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₃→∞ and k₁→0 while keeping k₁k₃ constant, in which case K₂→0 and K₁→∞ while K₂K₁ is constant (K₁ and K₂ were introduced as they appear in the 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₃. We used for positive cooperativity K₁=1000, K₂=0.001 and K₁=0.5, K₂=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.

Figure 11

Evolution of the product concentration in both positive and negative cooperativity.

In Figure 12, we plot the evolution of the enzyme concentration for both extreme positive cooperativity and negative cooperativity.

Figure 12

Evolution of the enzyme concentration in both positive 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 )=m a x{d i a m(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 . 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.

Table 1 Convergence results for the basic enzyme reaction

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.

Table 2 Convergence results for the 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.

Table 3 Convergence results for the 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, . 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.

Figure 13

Snapshots of the product concentration in the basic reaction at T=0.5 with three different time steps, shown in subfigures (a), (b) and (c). The time steps are given below each subfigure.

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 14

Snapshots of the product concentration in the suicide substrate reaction at T=4 with three different time steps, shown in subfigures (a), (b) and (c). The time steps are given below each subfigure.

Stability of the cooperative reaction

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.

Figure 15

Snapshots of the product concentration in the cooperative reaction at T=4 with three different time steps, shown in subfigures (a), (b) and (c). The time steps are given below each subfigure.

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–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.