A finite element method for growth in biological development

1.   Introduction

Ion channels are protein channels in the cell membrane that regulate the balance and transmission of ions inside and outside the cell. Ion channels are involved in critical physiological processes such as cellular potentiation, signaling, cellular excitability, and metabolic regulation by selectively allowing specific types of ions to pass through the cell membrane in order to maintain the difference between the ionic concentrations inside and outside the cell. In this way, ion channels control a wide range of biological functions of living organisms ^[1,2,3,4]. Currently, the two most relevant properties of ion channels, permeation and selectivity, are characterized by the I-V relations (current-voltage relations) for each ionic species ^[1,5]. The ionic transport through the ion channel is governed by fundamental physical laws of electrodiffusion. The macroscopic properties of ionic flows also depend on external driving forces, namely boundary concentrations and boundary potentials ^[6].

There are two related major topics in the study of ion channel problems: structures of ion channels and ionic flow properties. The physical structure of ion channels is defined by the channel shape and the spatial distribution of permanent and polarization charge. Usually, the shape of a typical ion channel is approximated as a cylindrical-like domain with a varying cross-section area. Within a large class of ion channels, amino acid side chains are distributed mainly over a "short" and "narrow" portion of the channel ^[7], with acidic side chains contributing permanent negative charges and basic side chains contributing permanent positive charges. For an open channel with a given structure, the main interest is to understand its electrodiffusion property. One major challenge in the study of ionic flow properties is the nonlinear interaction among different physical parameters involved in the diffusion process of charged particles. On the other hand, all present experimental measurements about ionic flow are of the input-output type ^[1]; that is, the internal dynamics within the channel cannot be measured with current technology. Thus, it is extremely hard to extract coherent properties or to formulate specific characteristic quantities from experimental measurements. Therefore, a suitable mathematical analysis with a physically sound model plays critical and unique roles for a possible comprehensive understanding of ion channel problems.

Poisson-Nernst-Planck (PNP) systems are basic primitive models for ionic motion through ion channels, in which the channel is considered as a line segment with a position-dependent conductivity related to the cross-sectional area. Recently, there have been some successes in mathematical analysis of PNP models for ionic flows through membrane channels. Those ^{[6,7,8,9,10,11,12,13,14,15,16,17]} that were studied in the framework of geometric singular perturbation analysis implied certain interesting threshold phenomena of ionic flows for relatively simple setups ^[7,18].

1.1.   The PNP model for ionic flows

Considering the structural characteristics of the channel, a basic continuum model for ionic flows is the PNP system, which treats the aqueous medium as a dielectric continuum (see ^{[19,20,21,22,23,24,25,26]}). PNP systems, under some further conditions, can be derived as reduced models from molecular dynamic models, Boltzmann equations, and variational principles ^{[27,28,29,30,31]}. The PNP type models have different levels of resolutions: The classical PNP system is the simplest one, which treats ions as point-charges and ignores ion-to-ion interactions ^{[32,33,34,35,36,37,38,39,40]}. The PNP-HS model takes into consideration volume exclusion by treating ions as hard-spheres with the charges located at their centers (see ^[8,12,41,42] and references therein). More sophisticated models, such as coupling PNP and Navier-Stokes equations for aqueous solutions (see, e.g., ^{[43,44,45,46,47,48,49]}), and a PNP-Cahn-Hilliard model which includes steric effects ^[50,51] have been developed. The sophisticated models are able to model the physical problem more precisely, however it is very difficult to analytically and even numerically study their dynamics. Here, we incorporate key features while keeping the equations amenable to both analysis and numerical simulations of ionic flows.

Since the channel is narrow, and one can effectively view it as a one-dimensional channel $[0, l]$, where $l$ is the length of the channel together with the baths linked by the channel. A quasi-one-dimensional steady-state PNP model for a mixture of $n$ ion species though a single channel reads ^[52]

where $e$ is the elementary charge, $k_B$ is the Boltzmann constant, $T$ is the absolute temperature, $\Phi(X)$ is the electric potential, $Q(X)$ is the permanent charge density of the channel, $\varepsilon_0(X)$ is the local dielectric coefficient, $\varepsilon_r(X)$ is the relative dielectric coefficient, and $A(X)$ represents the area of the cross-section over the point $X$ and, for the $j$th ion species, $C_j(X)$ is the concentration, $z_j$ is the valence, $\mu_j$ is the electrochemical potential, ${\cal J}_j$ is the flux density, and $D_j(X)$ is the diffusion coefficient.

For system (1.1), the following boundary conditions are posed ^[9] for $k = 1, 2, \cdots, n$

We point out that, for our analysis in current work, we are only interested in solutions of systems (1.1) and (1.2) with positive concentrations. For both time evolution of the three-dimensional PNP system and the quasi-one-dimensional PNP system with two ion species, the positivity-preserving properties for the ion concentration $C_k$ has been rigorously discussed in ^[53] (see Proposition $2.1$ for details). For the steady-state quasi-one-dimensional classical PNP system with $n\geq 3$ different types of ion species, the authors in ^[15] show that there is a unique solution of the system with $C_k > 0$. At the PDE level, we do not have a result when the system involves more than two different types of ion species. However, for the set-up in this work, it is not difficult to verify from ^[15] that the concentration $C_k$ is positive for $\nu$ sufficiently small. Interested readers can also refer to ^[51] for more discussions of the positivity-preserving properties of the concentrations from the point of view of a numerical approach. We also remark that our analysis is for the steady-state PNP model, and so the method in the current work cannot be directly applied to time-dependent PNP models and higher dimensional PNP models, but it provides important insights into related problems and could be a starting point for a more complicated model.

For a solution to the boundary value problem (BVP) (1.1) and (1.2), the current ${\cal I}$ is

For fixed boundary concentrations, the ${\cal J}_k$'s depend only on ${\cal V}$, and Eq (1.3) defines the so-called I-V relations.

1.2.   Bikerman's local hard-sphere model

The electrochemical potential $\mu_i(X)$ for the $i$th ion species consists of the ideal component $\mu_i^{id}(X)$ and the excess component $\mu_i^{ex}(X)$: $\mu_i(X) = \mu_i^{id}(X)+\mu_i^{ex}(X),$ where

where $C_0$ is some characteristic number density. The classical PNP system takes into consideration of the ideal component $\mu_i^{id}(X)$ only, which reflects the collision between ion particles and water (medium) molecules. Ion species in the classical PNP system are treated as point charges, and ion-to-ion interactions are ignored. The classical PNP system is reasonable in dilute cases. Meanwhile, this is a major weakness of the classical PNP model since many important properties of ionic flows do depend on finite ion sizes, particularly for those cations with the same valences but different finite ion sizes, such as Na$^+$ (sodium) and K$^+$ (potassium). This is closely related to the selectivity property of ion channels. To study the effects on ionic flows from ion sizes, one needs to consider excess components in the electrochemical potential, and one choice is to include hard-sphere (HS) potentials. PNP models with ion size effects have been investigated computationally with great successes (see ^{[24,28,30,44]}, etc.), and have been mathematically analyzed (see, for example, ^{[12,41,42,54,55]}).

In this work, we take Bikerman's local hard-sphere (LHS) model ^[56] to approximate $\mu_i^{ex}(X)$

where $\nu_j$ represents the volume of the $j$th ion species. Particularly, we take $\nu_n = \nu$ so that $\lambda_n = 1$ in our following discussion.

We point out that since $C_i$ is the number density of $i$th ion species, then, ${\sum_{j = 1}^nv_jC_j < 1}$. In this sense, Bikerman's LHS takes into consideration nonzero ion sizes.

1.3.   Further assumptions

For definiteness, we will take the following settings:

$({\rm{A1}}).$ Consider three ion species ($n = 3$) with $z_1 = z_2 = -z_3 = 1$.

$({\rm{A2}}).$ $Q(X) = 0$ over the whole interval $[0, 1]$.

$({\rm{A3}}).$ $\mu_i(X)$ consists of both the ideal component $\mu_i^{id}(X)$ and the LHS potential $\mu_i^{Bik}(X)$ in (1.5).

$({\rm{A4}}).$ $\varepsilon_r(X) = \varepsilon_r$ and ${\cal D}_i(X) = {\cal D}_i$.

Under the assumptions (A1)–(A4), system (1.1) reads

The boundary conditions become, for $i = 1, 2, 3,$

We first make a dimensionless rescaling following ^[57]. Set $C_0 = \max\{{\mathcal L}_i, {\mathcal R}_i:i = 1, 2\}$ and let

Correspondingly, for $k = 1, 2, 3$, the BVP (1.6) and (1.7) becomes

with the boundary conditions

We take $h(x) = 1$ over $[0, 1]$ in our analysis since eventually it appears in the form of $\int_0^1\frac{1}{h(s)}ds$, a positive constant, which will not affect the behavior of ionic flows.

2.   Materials and methods

We take advantage of the work done in ^[8], and further apply regular perturbation analysis to study the effects on the I-V relations from finite ion sizes under further restrictions. To be specific, we expand the I-V relation $I(V, \nu)$ along $\nu = 0$ for fixed boundary concentrations and obtain

where, under our assumptions,

Here, from ^[8],

where

and

3.   Results

Our main focus in this section is to examine the finite ion size effects on the I-V relations. To be specific, we are interested in the leading term $I_1(V)$, which contains the finite ion size effects. More precisely, with respect to the membrane potential $V$, we will characterize

$({\rm{i}})$ the monotonicity of $I_1(V)$;

$({\rm{ii}})$ the sign of $I_1(V)$;

$({\rm{iii}})$ the magnitude of $I(V)$ equivalent to $I_0(V)I_1(V)$.

Our analysis is further based on the following assumptions

where $L_3 = L_1+L_2$ and $R_3 = R_1+R_2$ are the so-called electroneutrality boundary conditions.

To get started, we rewrite $J_{k1}$ for $k = 1, 2, 3$, under assumption (3.1), as follows:

where

It then follows that

For $J_{k0}, \ I_0, \ J_{k1},$ and $I_1$, we observe the following:

Lemma 3.1. The following scaling law holds.

$\rm{(i)}$ $J_{k0}$ and $I_0$ scale linearly in the boundary concentrations; that is, for any real number $s > 0$,

and

$\rm{(ii)}\ J_{k1}$ and $I_1$ scale quadratically in the boundary concentrations; that is, for any $s > 0$,

and

We remark that the scaling law actually provides an efficient way to adjust the effects on ionic flows from the finite ion sizes by controlling the boundary concentrations.

3.1.   The monotonicity of $I_1(V)$

We study the monotonicity of the leading term $I_1(V)$ in (3.3) under the assumption (3.1), which is critical for one to further consider the sign of $I_1(V)$.

To get started, we establish the following result.

Lemma 3.2. Assume $L_3 < R_3$. One has $q_2 > 0$.

Proof. We rewrite $q_2$ as

Note that $\ln A(1) = \ln\frac{R_3}{L_3} > 0.$ One has that the sign of $q_2$ is the same as that of $(L_3+R_3)\ln A(1)+2(L_3-R_3).$ For convenience, we define

where

Upon introducing $x = \frac{L_3}{R_3}$, $\bar{q}_2(L_3, R_3)$ can be written as

Direct calculation gives $\bar{q}_2'(x) = -\ln x-\frac{1}{x}+1$ and $\bar{q}_2''(x) = \frac{1-x}{x^2} > 0$ for $0 < x < 1$, which implies that $\bar{q}'_2(x)$ is increasing in $x$. Note also that $\bar{q}_2'(1) = 0.$ One has $\bar{q}_2'(x) < 0$ for $0 < x < 1$, which indicates that $\bar{q}_2(x)$ is decreasing for $0 < x < 1$. Together with $\bar{q}_2(1) = 0$, one has $\bar{q}_2(x) > 0$ for $0 < x < 1$. Hence, $\hat{q}_2(x) > 0$ for $0 < x < 1$. Therefore, $q_2 > 0$ for $L_3 < R_3$.

We have the following result:

Theorem 3.3. Under the assumption (3.1), for $V > 0$, $I_1(V)$ is increasing in the potential $V$.

Proof. Under assumption (3.1), from (3.2), direct calculations yield

Note that, for $V > 0$,

Together with Lemma 3.2, one has $J_{11}'(V) > 0, \ J_{21}'(V) > 0,$ and $J_{31}'(V) < 0$. By the definition of $I_1(V) = D_1J_{11}(V)+D_2J_{21}(V)-D_3J_{31}(V)$, we have $I_1'(V) > 0$, and thus $I_1(V)$ is increasing in the potential $V$.

3.2.   Signs of $I_1(V)$

Based on the analysis in Section 3.1, we now study the sign of the leading term $I_1(V)$, which provides important information on the finite ion size effects of the I-V relations.

Lemma 3.4. Under assumption (3.1), one has

$\rm{(i)}$ If $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2) > 0$, then, $I_1 > 0$;

$\rm{(ii)}$ If $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2) < 0$, then, there exists a unique zero $V_{1s}^{zc}$ of $I_1(V) = 0$, such that $I_1 < 0$ for $(0, V_{1s}^{zc})$ and $I_1 > 0$ for $(V_{1s}^{zc}, +\infty)$.

Proof. Evaluating $I_1(V)$ in (3.3) at $V = 0$, one has

where

It is easy to verify that, under the assumption $\frac{L_1}{L_2} = \frac{R_1}{R_2}$, $K > 0$. Hence, the sign of $I_1(0)$ is determined by the factor $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2).$ More precisely,

● if $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2) > 0$, then $I_1(0) > 0$. Together with the fact that $I_1(V)$ is increasing, one has $I_1(V) > 0$ for all $V > 0$;

● if $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2) < 0$, then $I_1(V) = 0$ has a unique zero, say $V_{1s}^{zc}$ such that $I_1(V) < 0$ for $V\in(0, V_{1s}^{zc})$ while $I_1(V) > 0$ for $V\in(V_{1s}^{zc}, \infty)$ since again $I_1(V)$ is increasing in the potential $V$.

This completes the proof.

One of our main results then follows directly.

Theorem 3.5. Under assumption (3.1) and $V > 0$, one has, for small $\nu > 0$,

$\rm{(i)}$ if $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2) > 0$, then the finite ion size always enhances the current $I(V)$;

$\rm{(ii)}$ if $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2) < 0$, then the finite ion size enhances the current $I(V)$ if $V\in(V_1^{zc}, \infty),$ while it reduces the current $I(V)$ if $V\in (0, V_1^{zc})$.

3.3.   Effects on the magnitude of $I(V)$ from finite ion sizes

We consider the finite ion size effects on the magnitude of the I-V relation $I(V)$, that is $|I(V)|$, which is equivalent to analyzing $I_0(V)I_1(V)$. For the zeroth order term $I_0(V)$, one has

Lemma 3.6. Under assumption (3.1), one has $I_0(V)$ is increasing in the potential $V$; furthermore,

$\rm{(i)}$ If $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2) > 0$, then $I_0(V) > 0$;

$\rm{(ii)}$If $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2) < 0$, then there exists a unique zero $V_{0s}^{zc}$ of $I_0(V) = 0$, such that $I_0 < 0$ for $(0, V_{0s}^{zc})$ and $I_0 > 0$ for $(V_{0s}^{zc}, +\infty)$.

Proof. To study the monotonicity of $I_0(V)$, we will focus on the monotonicity of the individual fluxes $J_{k0}$ for $k = 1, 2, 3$ first. From (2.2), direct calculations yield

where, with $x = V$, $a = L_3/{R_3} > 0,$ and $b = L_1/{R_1} > 0$,

Since $g_0'(x) = (a-e^{-x})(b+e^{-x})e^x$, $g_0(x; a, b)$ has only one critical point $x^* = -\ln a$. It follows from $g_0(x^*; a, b) = 0$ and $\lim_{x\to\infty}g_0(x; a, b) = \infty$ that $g_0(x; a, b) > 0$ for $x\neq x^*$. Note that the sign of $\frac{d J_{10}}{dV}$ is the same as that of $g_0(x; a, b)$. One has that $J_{10}$ is increasing in the potential $V$. Similarly, one can show that $J_{20}(V)$ is also increasing in the potential $V$. As for $J_{30}(V)$, one has $\frac{dJ_{30}}{dV}(V; 0) = -zD_3f_0(L_3, R_3) < 0$ for all $V$, from which $J_{30}(V)$ is decreasing in the potential $V$. It then directly follows from the definition of $I_0 = D_1J_{10}+D_2J_{20}-D_3J_{30}$ that the zeroth order term $I_0(V)$ is increasing in the potential $V$. Also, from (2.2), we have

Note that $I_0(0) = (D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2).$ We have

● If $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2) > 0$, then $I_0(V) > 0$;

● If $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2) < 0$, then there exists a unique zero $V_0^{zc}$ of $I_0(V) = 0$, such that $I_0 < 0$ for $(0, V_0^{zc})$ and $I_0 > 0$ for $(V_0^{zc}, +\infty).$

This completes the proof. It follows that:

Theorem 3.7. Suppose $\frac{L_1}{L_2} = \frac{R_1}{R_2}$. For $V > 0$, one has

$\rm{(i)}$ For $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2) > 0$, one has $I_0(V)I_1(V) > 0$;

$\rm{(ii)}$ For $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2) < 0$, one has

$\rm{(a)}$ If $V_0^{zc} < V_1^{zc}$, then

$\rm{(a1)}$ $I_0I_1 < 0$ if $V\in(V_{0}^{zc}, V_1^{zc})$;

$\rm{(a2)}$ $I_0I_1 > 0$ if $V\in(0, V_{0}^{zc})\cup(V_1^{zc}, +\infty)$.

That is, for $V\in(V_{0}^{zc}, V_1^{zc})$, ion sizes $\nu$ reduces $|I(\nu)|$, and for $V\in(0, V_{0}^{zc})\cup(V_1^{zc}, +\infty)$, ion sizes $\nu$ strengthens $|I(\nu)|$.

$\rm (b)$ If $V_1^{zc}<V_{0}^{zc}$, then

$\rm{(b1)}$ $I_0I_1<0$ if $V\in(V_1^{zc}, V_{0}^{zc})$;

$\rm{(b2)}$ $I_0I_1>0$ if $V\in(0, V_1^{zc})\cup(V_{0}^{zc}, +\infty)$.

That is, for $V\in(V_1^{zc}, V_{0}^{zc})$, ion sizes $\nu$ reduces $|I(\nu)|$, and for $V\in(0, V_1^{zc})\cup(V_{0}^{zc}, +\infty)$, ion sizes $\nu$ strengthens $|I(\nu)|$.

To end this section, we demonstrate that the ionic flow properties of interest are very sensitive to the sign of $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2)$, the interplay between the boundary concentrations and the diffusion coefficients. For any experiment to be conducted, once the ion species are fixed, the values for the diffusion coefficients $D_1, D_2,$ and $D_3$ will be fixed, and one can then adjust/control the boundary concentrations to observe different dynamics of ionic flows through membrane channels, which further depends on the boundary potential.

4.   Numerical simulations

In this part, numerical simulations are performed to provide more intuitive illustrations of some analytical results. To be specific, we numerically identify the critical potentials $V_{0}^{zc}$ and $V_1^{zc}$, which characterize the effects on ionic flows from finite ion sizes. To further illustrate the finite ion size effects on ionic flows, we also numerically obtain the zeroth order I-V relations (corresponding to the classical PNP system which ignores the finite ion sizes) and the I-V relations with small positive $\nu$, respectively. Corresponding critical potentials for each setup are identified, from which one is able to observe the effects on ionic flows from the finite ion sizes.

We will conduct two experiments:

$({\rm{I}})$ For the first one, we perform numerical simulations to the system under the assumption $\frac{L_1}{L_2} = \frac{R_1}{R_2}$, from which one can better understand our analytical results obtained in Section 3 (see Figure 1);

$({\rm{II}})$ For the second one, we perform numerical simulations to the system without the restriction $\frac{L_1}{L_2} = \frac{R_1}{R_2}$, which provides more general insights into the study of ionic flows properties of interest (see Figure 2).

To get started, we rewrite system (1.9) and (1.10) as a system of first-order ordinary differential equations. Upon introducing $u = \varepsilon\frac{d}{dx}\phi$, one has

where

The boundary conditions are

In our simulations of system (4.1) and (4.2), we take $\varepsilon = 0.08, \ \lambda_1 = 1.32, \ \lambda_2 = 1.85$ and $\nu = 0.01$ for both experiments, while taking $L_1 = 12, \ L_2 = 8, \ R_1 = 24, \ R_2 = 16$ for the first one with $\frac{L_1}{L_2} = \frac{R_1}{R_2}$, and $L_1 = 12, \ L_2 = 8, \ R_1 = 20, \ R_2 = 12$ for the second one. For each experiment, we take $D_1 = 1.032, \ D_2 = 1.334, \ D_3 = 1.45$ and $D_1 = 1.032, \ D_2 = 1.334, \ D_3 = 0.98$, respectively.

Our numerical simulations show that

$({\rm{A}})$ Under the restriction $\frac{L_1}{L_2} = \frac{R_1}{R_2}$, both the zeroth order term $I_0(V)$ and the leading term $I_1(V)$ that contains the finite ion size effects, are increasing in the boundary potential $V$ (see Figure 1). For $V > 0$,

$({\rm{A1}})$ if the quantity $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2)$ is positive, then both $I_0(V)$ and $I_1(V)$ are positive, and it follows that $I_0(V)I_1(V) > 0$ for $V > 0$;

$({\rm{A2}})$ if the quantity $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2)$ is negative, then both $I_0(V) = 0$ and $I_1(V) = 0$ have unique zeros, $V_{0s}^{zc}$ and $V_{1s}^{zc}$. That is, $I_0(V_{0s}^{zc}) = 0$ and $I_1(V_{1s}^{zc}) = 0.$ For the setup in our numerical simulations, $V_{0s}^{zc} < V_{1s}^{zc}$, and hence, $I_0(V)I_1(V) > 0$ for $V\in(0, V_{0s}^{zc})\cup(V_{1s}^{zc}, \infty)$, and $I_0(V)I_1(V) < 0$ for $V\in(V_{0s}^{zc}, V_{1s}^{zc}).$

This is consistent with our analytical results stated in Theorem 3.3, Lemma 3.4, Lemma 3.6, and Theorem 3.7.

(B) Without the restriction $\frac{L_1}{L_2} = \frac{R_1}{R_2}$, for $V > 0$, similar numerical results are obtained (see Figure 2). This indicates the qualitative properties of ionic flows through membrane channels are stable.

To end this section, we demonstrate that the purpose of the numerical simulations is just to provide more intuitive illustration of our analytical results, and it is simply based on "bvp4c" from [MATLAB]. There is no numerical convergence analysis involved, which is also beyond the goal of this work. Those who are interested in the numerical approach to similar problems may refer to the following works ^[51,58].

5.   Conclusions

In this work, we examine the qualitative properties of ionic flows through membrane channels via PNP systems with two cations and finite ion sizes effects modeled through Bikerman's local hard-sphere potential. The main goal is to understand the problem from the mathematical point of view, which could provide some insights into related studies of ion channel problems. Viewing the ion sizes as small parameters, we expand the I-V relation along the characteristic ion size $\nu = 0$ as

where $I_1(V)$ is the leading term that contains the finite ion size effects, and is our main focus in the study. The sign and the monotonicity of $I_1(V)$ is rigorously analyzed, from which one is able to examine the finite ion size effects on ionic flows. Critical potentials $V_{0s}^{zc}$ and $V_{1s}^{zc}$ balancing ion size effects on the I-V relations are identified, and their critical roles played in the study of ionic flow properties are analyzed. Nonlinear interplays between system parameters are characterized particularly, the ionic flow properties are sensitive to the interplay between the diffusion coefficients and the boundary concentrations of the two cations through $(D_1-D_3)(L_1-R_1)+(D_2-D_3)(L_2-R_2),$ which provides some efficient ways to adjust/control the finite ion size effects on the I-V relations. Interesting scaling laws are observed, which is another way to adjust/control the effects on the I-V relations from the finite ion sizes. Numerical simulations are then performed to provide more intuitive illustrations of the analytical results and better understanding of the internal dynamics of ionic flows, while current technology cannot detect it in experiments.

We comment that the setup in this work is relatively simple (e.g., no permanent charges considered), and the study is the first step for analysis on more realistic and complicated models. On the other hand, the simple model studied in this work allows one to obtain a more explicit expression of the I-V relations in terms of system parameters, from which we are able to extract concrete information of the effects from small finite ion sizes. Furthermore, the detailed discussion may have great potential for the studies of synthetic ion channels.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Y. Wang thanks the Math Department at New Mexico Tech for the hospitality during his one-year visit to the department when the main part of this work is completed under the instruction of M. Zhang. The authors are partially supported by Simons Foundation (No. 628308) and NSF of China (No. 12172199).

Conflict of interest

The authors declare there is no conflict of interest. Mingji Zhang is a guest editor for Mathematical Biosciences and Engineering and was not involved in the editorial review or the decision to publish this article.