New insights into the effects of small permanent charge on ionic flows: A higher order analysis

1.   Introduction

Ion channels, proteins within cell membranes, are vital for cell communication, signal transformation, and coordinated activities ^[1,2]. They are defined by their shape and permanent charge. These channels typically resemble cylinders, with amino acid side chains concentrated in a short and narrow region. Acidic side chains contribute negative charges, while basic side chains add positive charges, determining the channel's permanent charge. Channel structures selectively permit specific ions and ease their diffusion across cell membranes ^[3,4,5,6,7].

Permeation and selectivity properties of ion channels are currently derived mainly from experimentally measured current-voltage (Ⅰ-Ⅴ) relations ^[2,5,8,9]. While individual fluxes convey more detailed information, they are costly and difficult to measure ^[6,10]. The Ⅰ-Ⅴ relation reflects the channel structure's response to ionic fluxes but is influenced by boundary conditions that drive ionic transport ^[11,12]. This multi-scale nature, with various physical parameters, grants the system great flexibility and diverse behaviors—a hallmark of natural devices ^[3]. However, this complexity also poses challenges in extracting meaningful insights from experimental data, especially given the limitations in observing internal dynamics.

The Poisson-Nernst-Planck (PNP) model stands out as one of the most commonly utilized mathematical frameworks for studying ion channels ^{[13,14,15,16,17,18,19,20,21]}. This model takes into account the interplay between structural characteristics and physical parameters, and researchers have extensively examined it using a geometric singular perturbation (GSP) approach ^[22,23]. Through the application of this approach, the PNP model can be simplified into an algebraic system referred to as the governing system. Analyzing this governing system unveils crucial properties of ion channels, providing valuable insights for informed design and optimization across various applications ^[24,25].

The effects of permanent charge on ionic flows have been investigated by several studies using the PNP model, with both analytical and numerical methods ^{[10,26,27,28]}. Liu et al. ^[10,27] examined the flux ratios and ion channel structures via PNP, and analyzed how they influence the fluxes, boundary concentrations, and electric potentials of the system. Other papers explored the reversal potential and permanent charge under unequal diffusion coefficients, and derived universal properties of the system ^{[29,30,31,32,33]} or numerically studied the permanent charge effects on flux ratios, revealing new phenomena and qualitative changes ^[34,35]. These studies enhance the understanding of the channel geometry and the role of permanent charge in ion channel dynamics.

Given the complex and multi-scale nature of the problem at hand, a comprehensive understanding of the interactions between permanent charges and boundary conditions cannot solely rely on analytical methods, especially across varying magnitudes of permanent charges. Therefore, this study adopts a combined approach of analytical insights and numerical methods to delve deeper into how permanent charges influence ionic flows in the presence of electric potentials. Leveraging previous analytical work, particularly from studies such as ^[16,36], which investigated flux ratios under different conditions, the focus narrows down to the flux ratio introduced in ^[27], examining its dependencies on permanent charges and electric potentials.

The methodology integrates rigorous analysis and numerical simulations to explore the impact of permanent charges on individual fluxes within fixed-shape open channels. Rigorous analysis uncovers essential biological properties and classifies distinct behaviors across various physical domains, especially in limiting or ideal scenarios. Conversely, numerical simulations extend these analytical findings into realistic parameter ranges, often unveiling additional phenomena. This approach allows a deeper dive into the intricacies of permanent charge effects, expanding upon previous analyses based on the PNP framework, which have revealed intriguing phenomena related to small permanent charges ^[16].

For the numerical simulations, Python along with the Numpy and Matplotlib libraries ^[37] are utilized. The computational code used in this study is publicly accessible through the author's GitHub repository at https://github.com/Hamid-Mofidi/PNP/tree/main/Q2contribution. This open access repository encourages collaboration and facilitates knowledge sharing among researchers interested in this field.

In this manuscript, we revisit the zeroth and first order solutions in permanent charge, as outlined in ^[16], to pave the way for higher-order analyses. The new findings and highlights of our studies in this manuscript are as follows:

(a) We derive analytical expressions for the intricate second order solutions, which are elaborated in Section 3.2 and form the backbone of our study.

(b) We study the effect of permanent charge and boundary concentrations on fluxes and Ⅰ-Ⅴ relations, estimate error and assess nonlinear effects for fluxes (explained in Section 4).

(c) In Section 5, we explore the higher order impact of permanent charge on flux ratios and analyze their dependencies on voltages and permanent charges.

In addition, the combination of our analytical (Section 3) and numerical investigations (Sections 4 and 5) shed light on both linear and quadratic solutions, providing novel insights and expanding our understanding beyond existing frameworks. These results serve as the foundation for further exploration and analysis in this study.

The paper follows this structure: Section 2 introduces the classical PNP model for ion channels and establishes a quasi-one-dimensional electro-diffusion model in Section 2.1, considering two types of ions with different charges and a simple distribution of permanent charge. Section 2.2 transforms the model into a dimensionless form for simplified analysis. Section 2.3 presents the governing system for the boundary value problem (BVP). In Section 3, the singular solutions in the presence of small permanent charge are analyzed, exploring higher-order effects. Sections 3.1 and 3.2 respectively delve into the zeroth, first, and second order solutions and their implications for system behavior. Notably, Section 3.2 introduces new analytical results for the second order solutions in $Q_0$. Section 4 provides computational outcomes for the first and second order solutions in $Q_0$ and numerically investigates the impact of permanent charge on fluxes and Ⅰ-Ⅴ relations, revealing the intricate interplay between permanent charge, boundary conditions, and channel geometry. Sections 4.1 and 4.2 respectively focus on the first and second order effects. In Section 5, we study the higher order effects of small positive permanent charges on flux ratios. Finally, Section 6 concludes the manuscript, summarizing the main results, discussing implications, and suggesting directions for future research.

2.   Classical PNP systems for ion channels: Setup and key results

PNP systems, essential for studying ionic flows, originate from molecular dynamic models ^[38], Boltzmann equations ^[39], and variational principles ^[40,41]. Advanced coupling with Navier–Stokes equations ^[42,43,44] and rigorous establishment of the Onsager reciprocal law ^[45] offer sophisticated insights, striking a balance between accuracy and analytical/computational challenges, supported by reviews and model comparisons ^[46,47].

Building upon this foundation, we further streamline PNP models, especially for ion channels with narrow cross-sections relative to lengths, resulting in quasi-one-dimensional models ^[48]. This reduction yields quasi-one-dimensional models ^[48], with rigorous justification provided in ^[49]. The streamlined approach addresses both accuracy and analytical/computational challenges.

This section provides a detailed exposition of our mathematical model for ionic flows, focusing on the essential setup and key results. Specifically, we explore a quasi-one-dimensional PNP model that characterizes ion transport within a confined channel featuring a permanent charge. To ensure clarity in our subsequent analysis, we introduce notation and assumptions consistently used throughout the paper. Moreover, we review relevant findings from previous literature, such as ^[14,18], serving as crucial foundations for our contributions outlined in the following sections.

Remark 2.1. The time-dependent PNP model has been discussed in ^[14]. Equation (2.1) in the following is selected based on two primary reasons: First, the one-dimensional system offers simplicity. Second, if the one-dimensional limiting system maintains structural stability, the dynamics of the three-dimensional system mirror those of the one-dimensional counterpart. Verification of structural stability follows a well-established framework, albeit nontrivial. Therefore, understanding the behavior of the steady-state in the limiting one-dimensional system serves as a pivotal step in this context.

2.1.   A quasi-one-dimensional PNP model

Our analysis is based on a quasi-one-dimensional PNP model first proposed in ^[48] and, for a special case, rigorously justified in ^[49]. For a mixture of $n$ ion species, a quasi-one-dimensional PNP model is

where $X\in [a_0, b_0]$ is the coordinate along the axis of the channel and baths of total length $b_0-a_0$, $A(X)$ is the area of cross-section of the channel over the longitudinal location $X$, $e_0$ is the elementary charge, $\varepsilon_0$ is the vacuum permittivity, $\varepsilon_r(X)$ is the relative dielectric coefficient, ${\mathcal{Q}}(X)$ is the permanent charge density, $k_B$ is the Boltzmann constant, $T$ is the absolute temperature, $\Phi$ is the electric potential, for the $k$th ion species, $C_k$ is the concentration, $z_k$ is the valence, ${\mathcal{D}}_k(X)$ is the diffusion coefficient, $\mu_k$ is the electrochemical potential, and ${\mathcal{J}}_k$ is the flux density.

Equipped with the system (2.1), a meaningful boundary condition for ionic flow through ion channels (see, ^[14] for reasoning) is, for $k = 1, 2, \cdots, n$,

In relation to typical experimental designs, the positions $X = a_0$ and $X = b_0$ are located in the baths separated by the channel and are locations for two electrodes that are applied to control or drive the ionic flow through the ion channel. An important measurement is the Ⅰ-Ⅴ (current-voltage) relation where, for fixed ${\mathcal{L}}_k$'s and ${\mathcal{R}}_k$'s, the current ${\mathcal{I}}$ depends on the transmembrane potential (voltage) ${\mathcal{V}}$ by ${\mathcal{I}} = \sum_{s = 1}^nz_s{\mathcal{J}}_s({\mathcal{V}}).$

Certainly, the relations of individual fluxes ${\mathcal{J}}_k$ with respect to ${\mathcal{V}}$ are more informative, but, measuring them experimentally is much more difficult ^[50]. Ideally, the experimental designs should not affect the intrinsic ionic flow properties so one would like to design the boundary conditions to meet the so-called electroneutrality $\sum_{s = 1}^nz_s{\mathcal{L}}_s = 0 = \sum_{s = 1}^nz_s{\mathcal{R}}_s.$ The reason for this is that, otherwise, there will be sharp boundary layers which cause significant changes (large gradients) of the electric potential and concentrations near the boundaries so that a measurement of these values has nontrivial uncertainties. One smart design to remedy this potential problem is the "four-electrode-design": two 'outer electrodes' in the baths far away from the ends of the ion channel to provide the driving force and two 'inner electrodes' in the baths near the ends of the ion channel to measure the electric potential and the concentrations as the "real" boundary conditions for the ionic flow. At the inner electrodes locations, the electroneutrality conditions are reasonably satisfied, and hence, the electric potential and concentrations vary slowly and a measurement of these values would be robust. The cross-sectional area $A(X)$ generally exhibits the characteristic of being significantly smaller for $X$ in the interval $(a_0, b_0)$ (representing the neck region of the channel) compared to $X$ outside the interval $[a_0, b_0]$.

2.2.   Dimensionless form of the quasi-one-dimensional PNP model

The following rescaling or its variations have been widely used for the convenience of mathematical analysis ^[51,52]. Let $C_0$ be a characteristic concentration of the ion solution. We now make a dimensionless rescaling of the variables in the system (2.1) as follows.

We assume $C_0$ is fixed but large so that the parameter $\varepsilon$ is small. Note that $\varepsilon = \lambda_D/(b_0-a_0)$, where $\lambda_D$ is the Debye screening length. In terms of the new variables, the BVP (2.1) and (2.2) becomes

with boundary conditions at $x = 0$ and $x = 1$

where $V: = \frac{e_0}{k_BT}{\mathcal{V}}, \quad L_k: = \frac{{\mathcal{L}}_k}{C_0}, \quad R_k: = \frac{{\mathcal{R}}_k}{C_0}.$ The permanent charge $Q(x)$ is

where $0 < a = \frac{A-a_0}{a_1-a_0} < b = \frac{B-a_0}{a_1-a_0} < 1.$ We will take the ideal component $\mu_k^{id}$ only for the electrochemical potential. In terms of the new variables, it becomes

The ideal component $\mu_k^{id}(x)$ contains contributions of ion particles as point charges and ignores the ion-to-ion interaction. PNP models including ideal components are referred to as classical PNP models. Recall that the critical assumption is that $\varepsilon$ is small. This assumption allows us to treat the BVP (2.4) with (2.5) as a singularly perturbed problem. A general framework for analyzing such singularly perturbed BVPs in PNP-type systems has been developed in prior works ^[14,18] for classical PNP systems and in ^[34,52,53] for PNP systems with finite ion sizes.

A method described in ^[14] (expanded upon in ^[18]) addresses the connection issue within classical PNP models by breaking down the system into two subsystems: the fast and slow systems under limiting conditions. Leveraging the specific structures of the PNP system allows for the integration of these subsystems, resulting in the creation of a singular orbit as an initial approximation. Aligning slow and fast orbits gives rise to a set of algebraic equations governing these singular orbits. This study takes a direct approach using regular perturbation theory to derive singular orbits for small magnitudes of $|Q_0|$, complementing the broader methodology detailed in previous literature.

We now recall the result from ^[14], upon which our work will be based. For $n = 2$ with $z_1 > 0 > z_2$, the authors applied the GSP theory to construct the singular orbit of the BVP (2.4) and (2.5). The BVP is then reduced to a connecting problem: finding an orbit from $B_0 = \{(V, u, L_1, L_2, J_1, J_2, 0):\; \text{arbitrary } u, J_1, J_2\},$ to $B_1 = \{(0, u, R_1, R_2, J_1, J_2, 1):\; \text{arbitrary } u, J_1, J_2\}.$

On each interval, a singular orbit typically consists of two singular layers and one regular layer:

(1) In view of the jumps of permanent charge $Q(x)$ at $x = a$ and $x = b$, the construction of singular orbits is split into three intervals $[0, a]$, $[a, b]$, $[b, 1]$, as depicted in Figure 1. To do so, one introduces (unknown) values of $(\phi, c_1, c_2)$ at $x = a$ and $x = b$:

These values then determine boundary conditions at $x = a$ and $x = b$ as $B_a = \{(\phi^a, u, c_1^a, c_2^a, J_1, J_2, a):\; \text{arbitrary } u, J_1, J_2\},$ and $B_b = \{(\phi^b, u, c_1^b, c_2^b, J_1, J_2, b): \; \text{arbitrary } u, J_1, J_2\}.$ Consequently, there are six unknowns $\phi^{a}$, $\phi^b$, $c_k^a$, and $c_k^b$ for $k = 1, 2$ that should be determined. On interval $[0, a]$, a singular orbit from $B_0$ to $B_a$ consists of two singular layers located at $x = 0$ and $x = a$, denoted as $\Gamma_0^l$ and $\Gamma_a^l$, and one regular layer $\Lambda_l$. Furthermore, with the preassigned values $\phi^a$, $c_1^a$, and $c_2^a$, the flux $J_k^l$ and $u_l(a)$ are uniquely determined so that $(\phi^a, u_l(a), c_1^a, c_2^a, J_1^l, J_2^l, a)\in B_a.$

(2) On interval $[a, b]$, a singular orbit from $B_a$ to $B_b$ consists of two singular layers located at $x = a$ and $x = b$, denoted as $\Gamma_a^r$ and $\Gamma_b^l$, and one regular layer $\Lambda_m$. Furthermore, with the preassigned values $(\phi^a, c_1^a, c_2^a)$ and $(\phi^b, c_1^b, c_2^b)$, the flux $J_k^m$, $u_m(a)$ and $u_m(b)$ are uniquely determined so that $(\phi^a, u_m(a), c_1^a, c_2^a, J_1^m, J_2^m, a)\in B_a\; \text{ and }\; (\phi^b, u_m(b), c_1^b, c_2^b, J_1^m, J_2^m, b)\in B_b.$

(3) On interval $[b, 1]$, a singular orbit from $B_b$ to $B_1$ consists of two singular layers are located at $x = b$ and $x = 1$, denoted as $\Gamma_b^r$ and $\Gamma_1^l$, and one regular layer $\Lambda_r$. Furthermore, with the preassigned values $\phi^b$, $c_1^b$, and $c_2^b$, the flux $J_k^r$ and $u_r(b)$ are uniquely determined so that $(\phi^b, u_r(b), c_1^b, c_2^b, J_1^r, J_2^r, b)\in B_b.$

2.3.   Governing system for the BVP

The matching conditions of the connecting problem in the previous section are

There are a total of six conditions, which are exactly the same number of unknowns preassigned in (2.8). Then, the singular connecting problem is reduced to the governing system (2.9) (see ^[14] for an explicit form of the governing system). More precisely,

where $y > 0$ is also unknown, and under electroneutrality boundary conditions $z_1L_1 = -z_2L_2 = L$ and $z_1R_1 = -z_2R_2 = R$,

Remark 2.2. In (2.10), the unknowns are: $\phi^{a}, \; \phi^{b}, \; c_1^{a}, \; c_2^{a}, \; c_1^{b}, \; c_2^{b}, \; J_1, \; \phi^{a, m}, \; \phi^{b, m}, \; y^*$, and $Q_0$, that is, there are eleven unknowns that match the total number of equations on (2.10).

Remark 2.3. In the upcoming sections, we will encounter lengthy terms in certain formulas. To simplify our notation, we introduce the following abbreviations for $k = 0, 1, 2$:

3.   Expanding singular solutions in the presence of small permanent charge

This section, and particularly Section 3.2, involves numerous intricate computations, undertaken with rigorous precision and validated through multiple verifications. However, to maintain readability, the detailed computations have been presented in compact form within the text. Interested readers are encouraged to meticulously examine each step and process to replicate the results accurately. Detailed computations pertaining to Section 3.1 can be found in the papers ^[16,17]. Additionally, for further clarification on Section 3.2, the authors are available upon request and can provide a comprehensive version of the paper to the journal if necessary.

Assuming that $|Q_0|$ is small, we expand all unknown quantities in the governing system (2.10) and (2.11) in $Q_0$, i.e., we write

where, $I_k$, for $k = 0, 1, 2$, were defined in (2.12).

Remark 3.1. To simplify matters, we made the assumption that all diffusion coefficients $D_k$ in (2.3) are equal. Therefore, we did not include them in our calculations in (3.1).

Remark 3.2. In the upcoming sections, as illustrated in (3.1), the subscripts 's' in $\phi_s^a, c_{ks}^a, J_{ks}$, etc., indicate the term's order when expanded with respect to $Q_0$. Here, 's' can represent values of $0, 1$, or 2, corresponding to the zeroth, first, or second-order term, respectively.

3.1.   Zeroth and first order solutions in $Q_0$ of (2.10) and (2.11).

The problem for the limiting case, where $Q_0 = 0$, has been addressed in ^[17] for $h(x) = 1$, and for a general $h(x)$, it can be resolved as demonstrated in ^[14] over the interval $[0, a]$. One can also derive the zeroth order solution directly by substituting Eq (3.1) into (2.10), expanding the identities in $Q_0$, and comparing the terms of like-powers in $Q_0$. Below, we outline the results for the zeroth and first order terms. The detailed proofs for these solutions can be referenced in ^[16]. These expressions are essential for the computational calculations discussed in Section 4.1, as well as for the computations related to second-order solutions in Section 4.2. Denote,

where $H(x)$ was defined in (2.11). Note that if $h(x)$ is uniform, then $H(x)$ represents the ratio of the length to the cross-sectional area of the portion of the channel over the interval from 0 to $x$ ^[35]. The origin of this quantity, $H(x)$, can be traced back to Ohm's law for the resistance of a uniform resistor. It is important to highlight that the parameters $a$ and $b$, along with the value $Q_0$, play pivotal roles in defining the shape and the permanent charge of the channel structure. For a more comprehensive understanding of the influences of $a$ and $b$ on the fluxes, refer to Section 4 in ^[16].

Proposition 3.1. The zeroth order solutions in $Q_0$ of (2.10) and (2.11), under electroneutrality boundary conditions $z_1L_1 = -z_2L_2 = L$ and $z_1R_1 = -z_2R_2 = R$ where one obtains $c_j^L = L_j, c_j^R = R_j, \phi^L = V, \phi^R = 0$, are given by

To compute the first-order terms in $Q_0$, we adopt the method introduced in ^[16], where we represent the intermediate variables in relation to the zeroth-order terms. The proof process is straightforward: by expanding the relevant identities in (2.11) with respect to $Q_0$, comparing the first-order terms in $Q_0$, and utilizing the results derived from Proposition (3.1), we can establish the desired relations.

Lemma 3.2. For the first order solution in $Q_0$ of (2.10) and (2.11), we obtain

By applying the same procedure as above to the remaining four identities in (2.10), and utilizing the results from Proposition (3.1) and Lemma (3.2), one can directly derive the first-order terms as follows.

Proposition 3.3. The first-order terms of the solution in $Q_0$ for the system (2.10) are as follows:

where, under electroneutrality boundary conditions $z_1L_1 = -z_2L_2 = L$ and $z_1R_1 = -z_2R_2 = R$, and in terms of $\alpha = \frac{H(a)}{H(1)}$ and $\beta = \frac{H(b)}{H(1)}$, the expressions for $A$, $B$ and $\lambda$ are

3.2.   Second order solutions in $Q_0$ of (2.10) and (2.11).

The results presented in this section extend the findings of the previous section, employing a consistent approach and methodologies with solutions exhibiting regularity concerning the permanent charge. To date, as far as we know, other papers have examined only up to the first-order terms in Eq (3.1), and the quadratic expression obtained in Section 3.2 is introduced for the first time in this work. Hence, all results in this section represent novel findings. Nevertheless, it is crucial to note that certain intricate calculations, owing to their extensive nature, are condensed for the sake of clarity in the presentation. For the second order solutions in terms of $Q_0$, we will first express the intermediate variables such as $\phi_{2}^{a, l}, \; c_{k2}^{a, l}$, etc. in terms of zeroth and first order terms and $\phi_{2}^a, \; c_{k2}^a$, etc.

Lemma 3.4. For the second-order solutions in terms of $Q_0$, we have

Proof. We present the derivations of the first two equations without showing the tedious computations, which mainly involve manipulating lengthy terms. The first step is to substitute (3.1) into the first equation in (2.10) and expand with respect to the parameter $Q_0$. Then, by applying a Taylor expansion for the function $e^{z_k(\phi^a - \phi^{a, m})}$ with respect to $Q_0$, we obtain the following expression for the second order terms:

Next, we substitute the expression for $c_1^{a, l}$ from $c_1^{a, l}$ from (2.11) into the third equation of (2.10) and expand the resulting equation up to third-order terms in $Q_0$, which gives us:

To obtain the desired result, we must carefully compute the expansions on both sides of (3.5) up to the third order and simplify the terms accordingly. to obtain the desired result.

Note that for small values of $Q_0$, we can make an approximation: $z_1c_1^a + z_2c_2^a \approx z_1c_{10}^a + z_2c_{20}^a = 0$, which implies that $-z_2c_2^a/z_1c_1^a \approx 1$. Moreover, in the proof, we applied the Maclaurin expansion of the natural logarithm, given by $\ln (x) = \ln (1 + (x-1)) = (x-1) - \frac{1}{2}(x-1)^2 +\cdots$. This expansion converges when $|x-1| < 1.$

Lemma 3.5. For the second-order intermediate variables in terms of $Q_0$, we establish

Proof. Starting from the second line of (2.11), we can derive the second order terms as follows:

By substituting $c_{11}^a$ from Proposition (3.3), we obtain the formula for $\phi_{2}^{a, l}$.

Moving on to the fourth line of (2.11), the second order terms can be expressed as:

Finally, from the sixth line of (2.11), we can determine $c_{12}^{a, m}$. Similar relations can be found for the other terms.

By following the previously outlined procedure for the last four identities in (2.10) and leveraging the results from Proposition (3.3), along with Lemmas (3.4) and (3.5), one can straightforwardly derive the following Lemma.

Lemma 3.6. Second order fluxes of the solution in $Q_0$ to the system 2.10 are given by

where,

and $T_0, T_1$ and $T_2$ were defined in (2.12).

Proof. Consider the expression $J_1$ in (2.10), and expand the terms with respect to $Q_0$ to get,

Therefore, the zeroth and first order terms in $Q_0$ of $J_1$ are,

The second-order term in $Q_0$ of $J_1$, with a careful computation, will be as follows,

By substituting $c_{11}^{a, l}$, $c_{12}^{a, l}$, and $\phi_{1}^{a, l}$ as per Lemmas 3.2 and 3.5, we can directly derive the formula for $J_{12}$. The expression for $J_{22}$ can be obtained in a similar manner.

Proposition 3.7. Second order intermediate concentration terms of the solution in $Q_0$ to the system 2.10 are given by

Proof. Initially, we start by adding up the expressions for $J_{12}$ and $J_{22}$ as outlined in the equations for $J_{12}$ and $J_{22}$ in Lemma 3.6, using careful simplification procedures. Afterward, we include $c_{22}^a$ and $c_{22}^b$ into the derived expression using the relevant expressions from Lemma 3.4. Through comprehensive computational analysis, we determine the expressions for $c_{12}^a$ and $c_{22}^a$.

In the process of determining the variable $y_2$, our initial step involves solving the equation for $c_{12}^{b}$ as presented in Lemma 3.6, specifically for $K_2$. Following this, we proceed to substitute the expressions for $K_1$ and $K_2$ and subsequently solve the equation for $y_2$, resulting in a simplified expression that provides the formula for $y_2$.

Remark 3.3. It is important to mention that because there are no explicit solutions for $\phi_2^a$ at this stage, the fluxes $J_{12}$ and $J_{22}$ in Lemma 3.6 cannot be expressed explicitly. Therefore, additional simplifications are required to compute the expression for $\phi_2^a$, as demonstrated in Proposition 3.8.

Utilizing the procedure outlined earlier, we shall extend our analysis to encompass the remaining four identities specified in Eq (2.10). With the foundational insights obtained from Proposition (3.1) and Lemma (3.2), we can proceed to systematically deduce the second order terms as delineated below.

Proposition 3.8. Under the electroneutrality boundary conditions, where $\phi^L = V, \phi^R = 0$, $z_lL_1 = -z_2L_2 = L$ and $z_lR_1 = -z_2R_2 = R$, the following results hold,

where,

and,

Furthermore, $z_1 c_1^L = z_1 L_1 = L, z_1 c_1^R = z_1 R_1 = R$ due to electroneutrality and $T_0, T_1$ were defined in (2.12).

Proof. Starting from the expressions for $J_{12}$ and $J_{22}$ derived in Lemma 3.6 and employing the relationships established in Lemma 3.5 and Proposition 3.7, and through meticulous computations, one can directly derive the second order terms for fluxes and electric potentials.

Remark 3.4. In Proposition (3.8), it is noteworthy that the following relationships hold:

wherein,

Remark 3.5. We emphasize once more the complex computations used to derive the second-order solutions in Section 3, although they are condensed for readability. Furthermore, Proposition 3.8 provides us with the necessary explicit expressions for the second-order solutions of the fluxes $J_{12}$ and $J_{22}$. However, given the complexity of these solutions, deriving additional analytical results to examine their impact on flux behavior would be very challenging. Hence, we turn to numerical investigations for further exploration in Section 4.

4.   Impact of permanent charge and boundary conditions on fluxes and Ⅰ-Ⅴ relations

In this section, we investigate how permanent charges and the boundary conditions impact the movement of individual fluxes and the current-voltage (Ⅰ-Ⅴ) relations. When the magnitude of $Q_0$ (a measure of permanent charge) is small, the flux ${J}_k$ for the $k$-th type of ion and the current ${I}$ can be represented as in (3.1). The quantities $J_{1k}$ and $J_{2k}$, where $k = 0, 1, 2$, capture the primary effects of permanent charges and channel shape on the flow of ions. We will analyze these quantities to understand their impact.

Remark 4.1. In the subsequent sections of this part, we conduct numerical simulations alongside an analysis of the equations in Section 3. The integration of numerical methods and analytical insights enhances our comprehension of the analytical findings. Specifically, the complexity of the quadratic solutions in Section 3.2 necessitates leveraging numerical observations to gain a deeper understanding of the second order solutions and their impact on the system. In our numerical simulations, we choose simplicity and specificity by setting $a = 1/3$, $b = 2/3$ in (2.6), and $h(x) = 1$. As a result, this yields $\alpha = 1/3$ and $\beta = 2/3$ in (3.2).

As highlighted in the Introduction section, our numerical methods are implemented using Python in conjunction with the Numpy and Matplotlib libraries. We create heatmaps, if possible, to visualize the signs of fluxes in various figures during our investigations, necessitating the identification of roots within the expressions. To accomplish this, we leverage combinations of Python functions, specifically np.where and np.isclose, for root finding purposes. Additionally, we also tried utilizing some other functions like the root function from the em scipy.optimize module, which is commonly used to solve systems of nonlinear equations. Due to the structured nature of our code, we prioritize the np.where and np.isclose functions. The np.isclose function, with its parameters rtol (relative tolerance), atol (absolute tolerance), and equal-nan (specifying 'Not a Number' (NaN) handling), generates a boolean array indicating element-wise equality within specified tolerances. Similarly, np.where(condition, [x, y]) selects elements from x or y based on a given condition, which proves useful when combined with np.isclose to locate indices satisfying the condition ^[37].

4.1.   Exploring first-order effects: Influence of permanent charges on fluxes

We begin by revisiting and simplifying specific findings from ^[16] and presenting numerical results for the first-order terms. Initially, we articulate Theorem 4.8 in ^[16], providing numerical insights, and subsequently expand on our findings based on further numerical investigations.

Suppose $B \neq 1$ where $B$ is defined as in (3.3). Let $V_q^1$ and $V_q^2$ be defined as follows:

Then the following cases arise:

(ⅰ) if $V_q^1 < 0 < V_q^2$, then, for $V > V_q^1$, a small positive $Q_0$ decreases $|J_1|$, and for $V < V_q^1$, it enhances $|J_1|$. Similarly, for $V > V_q^2$, a small positive $Q_0$ decreases $|J_2|$, and for $V < V_q^2$, it strengthens $|J_2|$; more precisely,

(ⅰ1) for $V \in (V_q^1, V_q^2)$, $J_{10}J_{11} < 0$ and $J_{20}J_{21} > 0$;

(ⅰ2) for $V < V_q^1$, $J_{10}J_{11} > 0$ and $J_{20}J_{21} > 0$;

(ⅰ3) for $V > V_q^2$, $J_{10}J_{11} < 0$ and $J_{20}J_{21} < 0$;

(ⅱ) if $V_q^1 > 0 > V_q^2$, then, for $V < V_q^1$, a small positive $Q_0$ decreases $|J_1|$, and for $V > V_q^1$, it enhances $|J_1|$. Similarly, for $V < V_q^2$, a small positive $Q_0$ decreases $|J_2|$, and for $V > V_q^2$, it strengthens $|J_2|$; more precisely,

(ⅱ1) for $V \in (V_q^2, V_q^1)$, $J_{10}J_{11} < 0$ and $J_{20}J_{21} > 0$;

(ⅱ2) for $V > V_q^1$, $J_{10}J_{11} > 0$ and $J_{20}J_{21} > 0$;

(ⅱ3) for $V < V_q^2$, $J_{10}J_{11} < 0$ and $J_{20}J_{21} < 0$.

The statements above are presented in a more streamlined form based on Theorem 4.8 in ^[16], which indicates that either case (ⅰ) or (ⅱ) may occur regardless of whether $L < R$ or $L > R$. It is crucial to highlight that the roots $V_q^1$ and $V_q^2$ in Eq (4.1) correspond to the roots of $J_{10}J_{11}$ and $J_{20}J_{21}$, respectively. This allows us to examine the effects of including linear terms $J_{11}$ or $J_{21}$. Now, let us extend the aforementioned findings based on the numerical observations depicted in Figure 2. It is also important to note that, with the simplifying assumptions in Remark 9, we have standardized the structure of channel geometry to focus on our primary objective, which is analyzing the relationships between the different orders of solutions in terms of $Q_0$ and boundary conditions, and their effects on fluxes and the Ⅰ-Ⅴ relation.

Now, we present a key observation derived from our numerical investigations, which sheds light on critical values and their impact on flux magnitudes under specific conditions regarding channel geometry. The findings from our analysis, summarized below, reveal significant insights into the behavior of fluxes $J_1$ and $J_2$ under varying voltage conditions and boundary concentrations:

Given boundary concentrations $L$ and $R$, and under certain conditions regarding channel geometry we have:

(ⅰ) There exist critical values $V_1^* < 0 < V_2^*$ such that:

(ⅰ1) For $V > V_1^*$, a small positive $Q_0$ reduces $|J_1|$, and for $V < V_1^*$, it amplifies $|J_1|$.

(ⅰ2) For $V > V_2^*$, a small positive $Q_0$ decreases $|J_2|$, and for $V < V_2^*$, it strengthens $|J_2|$.

(ⅱ) There exist critical points $V_2^* < 0 < V_1^*$ leading to:

(ⅱ1) For $V < V_1^*$, a small positive $Q_0$ reduces $|J_1|$, and for $V > V_1^*$, it amplifies $|J_1|$.

(ⅱ2) For $V < V_2^*$, a small positive $Q_0$ decreases $|J_2|$, and for $V > V_2^*$, it strengthens $|J_2|$.

It is worth nothing that the detailed cases can be articulated similarly to the detailed parts in Theorem 4.8 in ^[16]. Additionally, it is important to notice that $V_1^*$ and $V_2^*$ denote $V_q^1$ and $V_q^2$ respectively, although they are derived from numerical results.

Below is an observation derived from our numerical investigations. While the theoretical proof is not overly challenging, we opt not to consider it. The following findings are noted:

(a) When the boundary concentrations on the left and right ($L$ and $R$, respectively) are nearly identical ($L \approx R$), a small positive $Q_0$ decreases $|J_1|$ while increasing $|J_2|$. This behavior is exemplified by the blue region near $L = 1$ in Figure 2(A) and the red region near $L = 1$ in Figure 2(B).

(b) When the boundary concentrations are significantly different, the voltage ranges, for which a small positive $Q_0$ reduces $|J_1|$, become smaller, as do the voltage ranges for which $Q_0$ raises $|J_2|$. In other words, as the gap between $L$ and $R$ widens, the blue region in panel (A) and the red region in panel (B) in Figure 2 also reduce. In particular, for a fixed $R$, let $V_1^*$ and $V_2^*$ represent the critical points of $J_{10}J_{11}$ and $J_{20}J_{21}$ corresponding to $L_1$ and $L_2$, respectively, and let $\bar{V}_1^*$ and $\bar{V}_2^*$ denote the critical points of $J_{10}J_{11}$ and $J_{20}J_{21}$ corresponding to $\bar{L}_1$ and $\bar{L}_2$. Refer to Figure 2. This yields:

(b.1) if $R < L_1 < \bar{L}_1$, then $V_1^* < \bar{V}_1^* < 0$ and $0 < \bar{V}_2^* < V_2^*$.

(b.2) if $\bar{L}_1 < L_1 < R$, then $0 < \bar{V}_1^* < V_1^*$ and $V_2^* < \bar{V}_2^* < 0$.

Remark 4.2. One can check the author's GitHub repository for additional validation of the above observation as well as the upcoming figures in the subsequent sections. Similar situations can also be seen for a fixed $L$ there. The GitHub repository link is: https://github.com/Hamid-Mofidi/PNP/tree/main/Q2contribution.

Remark 4.3. In computational plots where it appears that $l$ and $r$ are equal, it is essential to note that they are very close but not precisely equal. This distinction is crucial based on the results.

The numerical results shown in Figure 2 validate the discussed scenarios, where the right boundary concentration $R$ is fixed at 1, while $L$ is varied between 0 and 2. Initially, this figure presents individual heatmaps illustrating the signs of $J_{10} J_{11}$ and $J_{20} J_{21}$ to clarify their respective flux changes. The red regions indicate areas where $J_{10}$ and $J_{11}$ in panel (A) (or, equivalently, $J_{20}$ and $J_{21}$ in panel (B)) share the same signs, while the blue regions denote areas where the signs are opposite. The color scheme can be interpreted as follows:

a. Red regions indicate areas where a (small) positive $Q_0$ reinforces $|J_1|$ or $|J_2|$.

b. Blue regions denote areas where a (small) positive $Q_0$ diminishes $|J_1|$ or $|J_2|$.

Thus far, the validation of our computational approach and analytical findings has been achieved by comparing them with Theorem 4.8 in ^[16], incorporating both zeroth and first-order terms. Our numerical analyses not only confirmed these findings but also provided additional insights into the first-order terms.

Remark 4.4. To validate our findings, we employed two approaches in our numerical analysis:

1) First, we computed $V_q^1$ and $V_q^2$ according to Theorem 4.8 in ^[16], as outlined at the beginning of this section. We then determined the signs on each interval.

2) Second, we numerically identified the roots $V_1^*$ and $V_2^*$ without explicitly computing $V_q^1$ and $V_q^2$.

We confirmed that the results are consistent for the first-order terms. The latter approach is particularly advantageous when incorporating the second order terms in the subsequent section, as obtaining roots analytically could be challenging. This will be our approach in the following section.

The intricate nature of the second order terms, specifically the fluxes $J_{12}$ and $J_{22}$ discussed in Section 3.2, necessitates numerical approaches to determine their roots. Therefore, we turn to Python, leveraging the Numpy and Matplotlib libraries, to perform calculations for zeroth, first, and the second order terms ^[37]. Additionally, numerical tools are employed to identify flux roots, facilitating the study of their signs across diverse regions.

The theoretical analysis of complex second order terms in equations provided in Proposition 3.8 is challenging. As a result, we use computational methods to explore how permanent charges affect ion movement and the membrane's electrical behavior, focusing on the current-voltage (Ⅰ-Ⅴ) relation. We analyze and compare these outcomes to scenarios without permanent charges, examining how these differences affect membrane performance. Then we study higher order contributions of permanent charges. Our numerical investigation delves into understanding the intricate interactions of permanent charges, shedding light on their influence on crucial electrical properties. Through this exploration, our aim is to advance our comprehension of the system's behavior and offer valuable insights to the academic community.

4.2.   Exploring the effects of second-order solutions and boundary conditions on fluxes

As of now, to the best of our knowledge, previous studies have only investigated up to the first-order terms in Eq (3.1) ^[16], and the quadratic expression obtained in Section 3.2 is introduced for the first time in this work. Incorporating these quadratic terms into the linear solutions will increase the accuracy of the solutions, although it was analytically challenging to derive them.

In this section, we delve into the implications of incorporating the $Q_0^2$ term into the expressions. Our primary focus is on investigating how the inclusion of $J_{k2}$ influences the linear estimation of $J_1$, represented as $J_{k0} + J_{k1}Q_0$. Additional comprehensive and noteworthy findings have been uncovered. Using heatmaps to examine the signs of $(J_{k0} + J_{k1} Q_0) J_{12}$ for $k = 1, 2$, along with the product $J_{k0} J_{k1} J_{k2}$, has revealed deeper insights, highlighting the unique impact of the $Q_0^2$ terms on the results.

In ^[15], the author demonstrates that the sign of the flux $J_k$ for $k = 1, 2$ remains unaffected by a permanent charge. In biological and chemical terms, the sign of the flux is dictated by the driving force (the gradient of electrochemical potential) rather than the structure (permanent charge $Q_0$) of the channel protein. However, the magnitude of $J_k$ is indeed influenced by $Q_0$. Referring back to 3.1, where for $k = 1, 2$, we have $J_k = J_{k0} + J_{k1}Q_0 + J_{k2}Q_0^2 + O(Q_0^3)$, we now focus on analyzing the impact of $J_{k2}$, the second order term of the flux, on the magnitude and behavior of $J_k$ using various approaches:

(1) Computing the product $(J_{k0} + J_{k1}Q_0) J_{k2}$ to observe the effects of $J_{k2}$ on the linear estimation of $J_k$:

1.ⅰ. If the product is positive, it suggests that the presence of $J_{k2}$ amplifies the effect of the linear term $J_{k0} + J_{k1}Q_0$. This implies that the magnitude of the linear estimation of $J_k$ will increase.

1.ⅱ. Conversely, if the product is negative, it implies that $J_{k2}$ dampens or counteracts the effect of the linear term $J_{k0} + J_{k1}Q_0$, resulting in a decrease in the magnitude of the linear estimation of $J_k$ (see Figure 4).

(2) Assessing the joint effects of $J_{k0}$, $J_{k1}$, and $J_{k2}$ through the product $J_{k0} J_{k1} J_{k2}$:

2.ⅰ. This product considers the interaction between all three coefficients $J_{k0}$, $J_{k1}$, and $J_{k2}$. If the product is positive, it indicates a reinforcement of the flux $J_k$ by $J_{k2}$, leading to an increase in the magnitude of $J_k$. Conversely, a negative product suggests a damping effect on $|J_k|$ due to the combined influence of $J_{k0}$, $J_{k1}$, and $J_{k2}$ (see Figure 5).

(3) Computing the product $J_{k0}(J_{k1} + J_{k2}Q_0)$ to observe the effects of small $Q_0$ on $J_k$:

3.ⅰ. If the product is positive (negative), it indicates that a small positive $Q_0$ will enhance (diminish) $|J_k|$. However, since this scenario is akin to computing the product $J_{k0}J_{k1}$ in the previous section where $Q_0$ is small, it can be disregarded.

In the subsequent discussion, we will review the aforementioned cases. However, prior to that, we utilize Figure 3 to showcase the transformative impacts of introducing the $Q_0^2$ term, transitioning the behavior from linear to nonlinear (quadratic).

Figure 3 illustrates that, transitioning from linear to quadratic, the sign of the flux $J_1$ never changes in the observed cases, as expected, while the magnitude of $J_1$ increases. However, this method has several limitations: its primary constraint is its representation of only specific cases, which may not be indicative of other scenarios. Furthermore, despite providing similar figures, descerning whether the quadratic term diminishes or amplifies the flux remains challenging. Another limitation is its inability to clearly illustrate how the flux behaves for very small values of $V$.

4.2.1.   Exploring the influence of $J_{k2}$ on the magnitude of linear estimation of the flux $J_k$

In this section, we delve into the intricate relationship between the second order flux component, $J_{k2}$, and the linear estimation of flux $J_k$, where $k = 1, 2$. Given that $J_k \approx J_{k0}+J_{k1}Q_0$, we examine the influence of the second order flux, $J_{k2}$, on the flux $J_k$ for $k = 1, 2$ by analyzing its effects on the linear estimate of $J_k$. To initiate our exploration, we construct plots of the product $(J_{k0} + J_{k1}Q_0) J_{k2}$, which effectively illustrates the impact of the second order flux $J_{k2}$ on the linear approximation of $J_k$. Referencing Figure 4 specifically enables a visual comprehension of these effects for both $k = 1$ and $k = 2$ scenarios.

We emphasize that due to the smallness of $Q_0$, the term $J_{k1}$ could be ignored, facilitating a direct calculation of $J_{k0}J_{k2}$ to evaluate the influence of $J_{k2}$ on the linear approximation magnitude of $J_k$, yielding the same results.

The following insights are obtained from numerical observations:

(a) When the left and right boundary concentrations ($L$ and $R$, respectively) are almost equal ($L \approx R$), a small positive $Q_0$ increases the magnitude of linear estimations for both $J_1$ and $J_2$. This trend is highlighted by the red region near $L = 1$ in Figure 4.

(b) In cases when the boundary concentrations are unequal, denoted as $L \neq R$, two critical voltages $V_1^*$ and $V_2^*$ emerge:

b.ⅰ. For voltages $V$ within the range $(V_1^*, V_2^*)$, a small positive $Q_0$ reduces the magnitude of linear estimations for both $J_1$ and $J_2$.

b.ⅱ. Conversely, for voltages $V$ outside the range $V_1^*$ to $V_2^*$, i.e., for $V < V_1^*$ or $V > V_2^*$, a small positive $Q_0$ increases the magnitude of linear estimations for both $J_1$ and $J_2$.

Furthermore, as the difference between $L$ and $R$ increases, the voltage ranges where a small positive $Q_0$ reduces the magnitude of linear estimations for both $J_1$ and $J_2$ also expand.

4.2.2.   Evaluating the combined impact of $J_{10}$, $J_{11}$, and $J_{12}$ through their product

In this part, we explore the combined impact of $J_{10}$, $J_{11}$, and $J_{12}$ by examining their product, $J_{10} J_{11} J_{12}$. This product captures the interplay among all three coefficients: $J_{k0}$, $J_{k1}$, and $J_{k2}$. A positive product signifies a reinforcement of the flux $J_{k}$, resulting in an amplified magnitude of $J_k$. Conversely, a negative product indicates a damping effect on $|J_k|$, reflecting the combined influence of $J_{k0}$, $J_{k1}$, and $J_{k2}$.

The observations can be summarized as follows:

(a) When the left and right boundary concentrations are in proximity (i.e., $L \approx R$), there exists a single critical value $V_0^* \approx 0$. Under this condition:

a.ⅰ. If $V > V_0^*$, then $J_{k0}J_{k1}J_{k2} < 0$, leading to a damping effect on $|J_{k}|$ for $k = 1, 2$.

a.ⅱ. If $V < V_0^*$, then $J_{k0}J_{k1}J_{k2} > 0$, resulting in an amplified effect on $|J_{k}|$ for $k = 1, 2$. (Refer to Figure 5 near $L = 1$.)

(b) Conversely, when the left and right boundary concentrations are not sufficiently close, there are two critical voltages: one is $V_0^* \approx 0$, and the other could be either $V_-^* < 0$ or $V_+^* > 0$, depending on the channel geometry and boundary concentration values. Under these conditions:

b.ⅰ. If $V_0^*$ and $V_-^*$ are the critical values, then for $V$ in the interval $(V_-^*, V_0^*)$, $J_{k0}J_{k1}J_{k2} > 0$, leading to an increasing $|J_k|$, and for $V$ outside this interval, $J_{k0}J_{k1}J_{k2} < 0$, resulting in a diminishing effect on $|J_k|$.

b.ⅱ. However, if $V_0^*$ and $V_+^*$ are the critical values, then for $V$ in the interval $(V_0^*, V_+^*)$, $J_{k0}J_{k1}J_{k2} < 0$, causing decreasing of $|J_k|$, and for $V$ outside this interval, $J_{k0}J_{k1}J_{k2} > 0$, leading to a strengthened effect on $|J_k|$.

Additionally, as the disparity between $L$ and $R$ widens, the two critical voltages converge (refer to Figure 5 far from $L = 1$).

4.3.   Estimating errors and assessing nonlinear effects for fluxes

In the context of Taylor expansions and polynomial approximations, error estimation serves to assess the impact of truncating the series at a finite order. Neglecting higher-order terms in the expansion introduces approximation errors, which can lead to deviations from the true function behavior. Therefore, quantifying these errors is essential for ensuring the validity of the approximation and understanding its limitations. We explore the calculation and analysis of approximation errors introduced by neglecting higher-order terms in the expansion. By examining the magnitude and significance of these errors, we gain insights into the accuracy of the approximations and the necessity of including additional terms in the expansion.

We present a detailed methodology for computing and analyzing approximation errors in the context of Taylor expansions. This involves calculating the error term introduced by neglecting the $J_{12}$ and $J_{22}$ terms in the Taylor expansion:

This error estimation evaluates the approximate value missed by fluxes when using the linear term. In other words, incorporating $J_{12} Q_0^2$ allows us to approach the exact value of flux $J_1$, while omitting it provides an estimation of the error. Similarly, this applies to flux $J_2$ and the term $J_{22} Q_0^2$. Since $Q_0$ is small, $Q_0^2$ and, consequently, $J_{12} Q_0^2$ are also small. Therefore, we primarily focus on its sign to determine if flux $J_1$ gains a small value (if $J_{12}$ is positive) or misses a small value (if $J_{12}$ is negative). Figure 6(A) demonstrates that for a fixed $R$, there exists a voltage $V_L$ for any $L$, such that for $V < V_L$, flux $J_1$ misses, and for $V > V_L$, it gains a small value when the approximation becomes more accurate by adding the nonlinear term $J_{12}$. Similarly, for $J_2$, Figure 6(B) indicates that for a fixed $R$, there exists a voltage $V_L$ (same as for $J_1$) for any $L$, such that for $V < V_L$, flux $J_2$ gains, and for $V > V_L$, it misses a small value when the approximation becomes more accurate by adding the nonlinear term $J_{22}$. It is important to note that $L$ and $R$ are independent here; thus, if one fixes $L$, and varies $R$, a similar discussion applies (refer to the figures available in the GitHub repository).

Nevertheless, error estimation in Taylor expansions is not without its challenges. One significant challenge lies in finding the exact solution against which to compare the approximated results. In many cases, the exact solution may be unknown or difficult to determine, leading to uncertainties in the accuracy of the error estimation. Additionally, the choice of expansion point in the Taylor series can significantly impact the magnitude and behavior of the error. Moreover, there exists a trade-off between accuracy and computational cost when considering the inclusion of higher-order terms in the expansion. While including more terms may improve the accuracy of the approximation, it also increases the computational complexity and resource requirements. Acknowledging these obstacles underscores the complexity and uncertainty inherent in error estimation methods and highlights the need for further research to address these limitations and enhance the reliability of numerical approximations.

4.4.   Impact of higher-order solutions on the Ⅰ-Ⅴ relation

This section explores the impact of permanent charge on the current-voltage (Ⅰ-Ⅴ) relationship. We investigate how the presence of permanent charge influences the electrical behavior of the system. By examining the influence of permanent charge on electrical behavior, a deeper understanding of the underlying mechanisms governing charge transport and device performance is sought.

In Section 4.4 of ^[16], an analysis of the impact of small permanent charges $Q_0$ on Ⅰ-Ⅴ relations was conducted, specifically focusing on the zeroth and first-order terms of $Q_0$. Here, we present some of their findings concerning equal diffusion coefficients and subsequently validate these results through our numerical investigations. Furthermore, we extend the inquiry to higher orders and explore the influence of the second order terms in permanent charge $Q_0$ on the Ⅰ-Ⅴ relation.

The authors of ^[16] demonstrated that the current $I_0$ remains unaffected by boundary concentrations, while $I_1$ does depend on them. Specifically, they established the existence of a $V_{rev}$ (which equals zero when the diffusion coefficients $D_i'$s are equal) where $I_0 < 0$ if $V < V_{rev}$ and $I_0 < 0$ if $V > V_{rev}$, irrespective of $L$ and $R$ values. This finding aligns with our numerical investigations depicted in Figure 7(A). It is worth noting that when $V_{rev} = 0$ (for equal diffusion coefficients), the righthand sides of Figure 7(A) are red, indicating $I_0 > 0$ for $V > V_{rev} = 0$, and the left-hand side is blue, indicating $I_0 < 0$ for $V < V_{rev} = 0$, regardless of $L$ and $R$ values.

Additionally, Theorem 4.14 in ^[16] asserts that $I_1$ does vary with $L$ and $R$; specifically, for equal diffusion coefficients, $I_1 > 0$ when $L < R$ and $I_0 < 0$ when $L > R$. This result is also consistent with our numerical investigations shown in Figure 7(B) where $R = 1$ is fixed. The lower part of the figure, representing $L < R$, is red, indicating $I_1 > 0$, while the upper part, where $L > R$, is blue, indicating $I_1 < 0$ in this scenario.

Remark 4.5. In ^[16], it is additionally shown how the current $I_1$ is influenced by the boundary $V$ (in addition to $L$ and $R$) for nonequal diffusion coefficients. Nevertheless, we chose not to delve into the numerical results for this scenario in order to maintain our primary focus on the higher-order terms involving the permanent charge $Q_0$ and its impact on the Ⅰ-Ⅴ relation.

We now expand upon the findings regarding the second order terms of $Q_0$ using numerical analysis, as depicted in Figure 7(C). These investigations reveal that, akin to $I_0$, the current $I_2$ remains unaffected by the boundary concentrations $L$ and $R$ but is contingent solely upon the boundary $V$, which equals zero when diffusion coefficients are equal. Consequently, we derive the following relationship:

Remark 4.6. It is noteworthy to mention that in Eq (4.2), we could potentially assert $V > V_{rev} = 0$ (and similarly for the converse case).

In the following discussion, we employ Figures 8 and 9 (for fixed boundary concentrations at $L = 0.008, R = 0.001$, while varying $Q_0$ and $V$) to illustrate the significant effects of introducing the $Q_0^2$ term, which leads to a transition in behavior from linear to nonlinear (quadratic). Several key observations from these figures are outlined below:

(a) Monotonicity in $V$ for small $Q_0$ and non-monotonic behavior in $Q_0$ for fixed $V$: Observing Figure 8, we note that for small $Q_0$, the current $\mathcal{I}$ shows a monotone (increasing) behavior with respect to $V$. This aligns with the expectation that the current is primarily influenced by $V$ when $Q_0$ is small. However, as depicted in Figure 9, monotonicity does not hold concerning $Q_0$ when $V$ is held constant.

(b) Bifurcations of $\mathcal{I} = 0$ (reversal potential): For values of $Q_0$ near $Q_0 = 0.1$, as depicted in Figure 8, the second-order solution in terms of $Q_0$ for the current, $\mathcal{I}$, appears to have three roots. This phenomenon was not anticipated or predicted in earlier studies such as those mentioned in ^[16,54].

(c) Figure 8 ($Q_0 = 0.1$) shows that the dashed green curve for the quadratic current solution lies below the linear solution between two red circles representing negative and positive voltages. This indicates that a small positive $Q_0$ may increase or decrease the magnitude of linear current estimations.

5.   Permanent charge effects on flux ratios

This section delves into the influence of a positive permanent charge on the fluxes of both cation and anion species. To quantify this influence, we introduce the flux ratio $\lambda_k(Q_0) = J_k(Q_0)/J_k(0)$, which compares the flux $J_k(Q_0)$ associated with a permanent charge $Q_0$ to the flux $J_k(0)$ with zero permanent charge, for a given ion species under specific boundary conditions and channel geometry. Note that in the following, we may use $\lambda_k(Q_0)$, $\lambda_k(Q_0, V)$, or $\lambda_k(Q_0, Q_0^2)$ for simplicity or to demonstrate the dependence of flux ratios on $V$ or $Q_0^2$, respectively.

For $n = 2$ with $z_1 = 1$ and $z_2 = -1$, detailed analysis of the impact of permanent charge described by Eq (2.9) on flux ratios has been conducted for both small and large $Q_0$ ^[16,27,36]. The flux ratio $\lambda_k(Q_0)$ serves as a metric for measuring the impact of the permanent charge $Q_0$: When $\lambda_k(Q_0) > 1$, the flux is augmented by $Q_0$, and when $\lambda_k(Q_0) < 1$, the flux is diminished by $Q_0$. An analysis of PNP models governing ionic flows reveals a universal principle regarding the effects of permanent charge^[27]:

Proposition 5.1. For a positive permanent charge $Q_0$, if $\lambda_1(Q_0)$ denotes the flux ratio for cation species and $\lambda_2(Q_0)$ signifies the flux ratio for anion species, then $\lambda_1(Q_0) < \lambda_2(Q_0)$ holds true irrespective of boundary conditions and channel geometry.

This proposition is precise in that, especially for a small positive $Q_0$, various scenarios emerge based on boundary conditions and channel geometry, such as $(i)$ $\lambda_1(Q_0) < 1 < \lambda_2(Q_0)$, $(ii)$ $1 < \lambda_1(Q_0) < \lambda_2(Q_0)$, and $(iii)$ $\lambda_1(Q_0) < \lambda_2(Q_0) < 1$. Each of these options captures unique details in how flux ratios change with a positive permanent charge.

In the preceding section, we generated heatmaps to visualize flux sign patterns across various figures in our investigations, necessitating the identification of roots within the expressions. However, due to the computational complexity involved, not all initial values are conductive to heatmap creation. The extensive computations introduce numerous tiny errors across a range of values, ultimately resulting in significant final errors in the output, thus reducing the effectiveness of heatmaps for representation. Consequently, we refrain from employing heatmaps in the subsequent section and instead focus on specific cases to derive results.

5.1.   Dependence of flux ratios on $Q_0$ for fixed $V$

For small $Q_0$, as $V$ increases, the change occurs from $1 < \lambda_1 < \lambda_2$ to $\lambda_1 < 1 < \lambda_2$ and to $\lambda_1 < \lambda_2 < 1$ or from $\lambda_1 < \lambda_2 < 1$ to $\lambda_1 < 1 < \lambda_2$ and to $1 < \lambda_1 < \lambda_2$ (refer to ^[54]).

(ⅰ) Non-monotonic behavior in $Q_0$ for fixed values of $V$: Examining Figure 10 for $V = -10$ (and similarly for $V = 30$), it becomes apparent that while the flux ratio $\lambda_2$ appears to increase (or decrease) based on linear estimates, the quadratic approximations reveal a different trend.

(ⅱ) Possible Pitchfork Bifurcations at $\lambda_k = 1$: In Figure 10 (for $V = 30$), the behavior of $\lambda_k(Q_0, Q_0^2)$ for various $Q_0$ values, influenced by second-order solutions in $Q_0$, exhibits non-monotonic trends and can cross the value of 1 twice. This implies the existence of two instances where $\lambda_k(Q_0, Q_0^2) = 1$, or there is a $Q_0^*$ for which $\frac{\partial \lambda_k}{\partial V}(Q_0^*, V) = 0$, indicating the possibility of bifurcations at $\lambda_k(Q_0^*, V) = 1$. As we will see in Section 5.2 (part ii), there is a $V^*$ such that $\frac{\partial \lambda_k}{\partial V}(Q_0, V^*) = 0$. Consequently, there is a chance that the corresponding value of $Q_0$ for which $\frac{\partial \lambda_k}{\partial V}(Q_0, V^*) = 0$ is precisely $Q_0^*$. This implies the existence of $(Q_0^*, V^*)$ such that $\frac{\partial \lambda_k}{\partial V}(Q_0^*, V^*) = \frac{\partial \lambda_k}{\partial Q_0}(Q_0^*, V^*) = 0$, potentially leading to a pitchfork bifurcation.

5.2.   Dependence of flux ratios on $V$ for fixed permanent charges

We now examine the dependence of ion fluxes $\lambda_1$ and $\lambda_2$ on $V$ for several fixed values of $Q_0$. In Figure 11, they are plotted as functions of $V\in (-50, 50)$ for $Q_0 = 0.001$, $0.005$, and $0.01$.

(ⅰ) Monotonicity in $V$ for small $Q_0$. From Figure 11, one observes that for very small $Q_0 = 0.001$, both $\lambda_1$ and $\lambda_2$ are monotone (decreasing) in $V$. This is consistent with the theoretical prediction made in ^[16] and also with the numerical observations in ^[54], and with the intuition that the flux ratios are dominated by the effects of $V$ when $Q_0$ is small.

(ⅱ) Bifurcations at $\lambda_k = 1$: In Figure 11, for various $Q_0$ values, $\lambda_1(Q_0, Q_0^2)$, influenced by second-order solutions in $Q_0$, exhibits a discontinuity near $\lambda_1 = 1$, resulting in three values where $\lambda_1 = 1$. Additionally, with increasing $Q_0$, fluxes show non-monotonic behavior in $V$, crossing $\lambda_1 = 1$ multiple times. These behaviors were not predicted by the analysis in ^[16], but non-monotonicity was discussed in ^[54] through numerical observations. Similar discussions apply to $\lambda_2 = 1$.

6.   Concluding remarks and future work

In this study, we presented a comprehensive exploration of ion channel dynamics, focusing on the intricate influence of permanent charges. Theoretical and numerical analyses have been combined to unveil the qualitative shifts in fluxes, flux ratios, and electric potentials at higher-order contributions of permanent charge. The investigation has delved into the subtle interplay between boundary conditions and channel geometry, elucidating the nuanced impact of permanent charges on ion channel behavior. Our findings contribute to the understanding of ion electro-diffusion, shedding light on the complex interactions that arise due to permanent charges. The systematic perturbation analysis, spanning zeroth, first, and second order solutions, has provided valuable insights into the behavior of the system under the influence of small permanent charges. As we conclude this study, avenues for further research emerge.

As indicated in Remark 4.1, the complexity of the quadratic solutions in Section 3.2 led us to utilize numerical observations to further explore the second-order solutions and their effects on the system. Integrating numerical observations with analytical insights improves our understanding of the analytical results, which we plan to continue studying in the future. Additionally, the application of advanced numerical techniques and simulations may offer a more detailed understanding of ion channel behavior in complex biological environments. Further investigations could also delve into the impact of permanent charges on specific ion channel types, allowing for a more targeted analysis of their behavior. Moreover, experimental validation and comparison with existing biological data would provide a bridge between theoretical insights and real-world observations, enhancing the practical relevance of our findings.

Exploring local hard-sphere PNP systems, which account for finite ion sizes, offers valuable insights into the dynamics of ionic channels by considering ion sizes $d$ ^[53]. However, the computations become more complex in this case. A fascinating aspect of this study involves investigating higher-order solutions concerning ion size $d$ and permanent charge $Q_0$, specifically deriving $Q_0^2, Q_0d,$ and $d^2$ solutions. We derived solutions involving $Q_0^2$ in this manuscript. The work presented in ^[25] delves into the higher-order effects of ion size and provides $d^2$ solutions. Additionally, the paper ^[53] examines PNP models with ion size and permanent charge, and to complete the puzzle, one must carefully derive $Q_0d$ terms from that paper. By assembling all these quadratic terms, a more accurate exploration of the higher-order impacts of ion size and permanent charge becomes possible.

Use of AI tools declaration

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

Acknowledgments

The author acknowledges the valuable comments and support provided by Weishi Liu from the University of Kansas, Bob Eisenberg from Rush University Medical Center, and Mingji Zhang from New Mexico Tech University. Additionally, the author expresses gratitude for the valuable comments and support received from the anonymous reviewers.

Conflict of interest

The author declares there is no conflict of interest.