Modeling and Simulation of Thermo-Fluid-Electrochemical Ion Flow in Biological Channels

In this article we address the study of ion charge transport in the biological channels separating the intra and extracellular regions of a cell. The focus of the investigation is devoted to including thermal driving forces in the well-known velocity-extended Poisson-Nernst-Planck (vPNP) electrodiffusion model. Two extensions of the vPNP system are proposed: the velocity-extended Thermo-Hydrodynamic model (vTHD) and the velocity-extended Electro-Thermal model (vET). Both formulations are based on the principles of conservation of mass, momentum and energy, and collapse into the vPNP model under thermodynamical equilibrium conditions. Upon introducing a suitable one-dimensional geometrical representation of the channel, we discuss appropriate boundary conditions that depend only on effectively accessible measurable quantities. Then, we describe the novel models, the solution map used to iteratively solve them, and the mixed-hybrid flux-conservative stabilized finite element scheme used to discretize the linearized equations. Finally, we successfully apply our computational algorithms to the simulation of two different realistic biological channels: 1) the Gramicidin-A channel considered in~\cite{JeromeBPJ}; and 2) the bipolar nanofluidic diode considered in~\cite{Siwy7}.


Introduction and Motivation
Ion channels are ubiquitous in the cells of the human body. They allow communication between the intra and extra-cellular sites and are responsible for the signaling pathways that regulate the activity of each part of the body system. A complete presentation of ionic channels is far beyond the scope of this article. For a treatment of the general properties of ionic channels and of their characterization in cellular biology and neuroscience we refer to the books [19] and [15]. In broad terms, ion channel are "biological gates" that can be activated (opened) or deactivated (closed) by the application of electrical, chemical, mechanical or thermal stimuli.
In the present article we mainly focus on the interaction among electrical, chemical and thermal driving forces, starting from the evidence that the sensory system in mammals is capable of discriminating thermal stimuli ranging from noxious cold (< 8 ∘ C) to noxious heat (> 52 ∘ C). Of particular relevance are the biological temperature sensors denoted thermo-transient receptor potential (thermo-TRP) cation channels. Thermo-TRP channels constitute a "superfamily" of receptors in sensory neurons that can be activated by noxious stimuli, resulting in the transmission to the spinal cord and brain of the information of local pain perception [16]. The fact that the thermal thresholds of many thermo-TRP channels can be modulated by extracellular mediators has recently led to the development of antagonists for noxious channel blockage in therapeutic uses as novel analgesics [44]. This is the case, for example, of heat-sensitive TRP channels exposed to capsaicin, a natural ingredient of spicy foods such as red hot chilli peppers and of cold-sensitive TRP channels exposed to men-thol. In [9] and [46], it is shown how the chemical agonists capsaicin and menthol function for both types of TRP channels as gating modifiers, shifting activation curves towards physiological membrane potentials.
Despite a significant amount of experimental data, the mechanisms underlying the marked temperature sensitivity of channel gating are still largely unknown. For this reason, the use of mathematical tools is increasingly becoming popular to support biophysical conjectures and suggest novel theories for data interpretation. Two main kinds of approaches are presently utilized for channel modeling and simulation, one is called particle approach the other continuum approach. The particle approach is based on the assumption that ions are charged particles characterized by a given radius and mass. Each particle is in motion within the electrolyte fluid permeating the channel and interacts with the other particles according to Coulombian electrostatic forces directed along the radius vector connecting particle i with particle j. The corresponding simulation approach is based on the use of the Monte Carlo computational method in which the temporal dynamics of the ion population is studied by solving the motion equations self-consistently with the Poisson equation that provides the updated electric field at each position of the medium. For details on this procedure we refer to [6] and bibliography cited therein. The continuum approach is based on the assumption that each ion species constitutes an homogenized charged fluid. Each ion fluid flows throughout the electrolyte fluid according to the general laws of hydrodynamics to which the Poisson equation must be added to update the electric field at each position of the medium. The corresponding simulation approach is based on the use of the numerical methodologies that are commonly employed to treat the equations for the dynamics of a fluid, typically the finite element method or the finite volume method with some specific stabilization tools for handling the presence of strong convective terms. For details on this procedure we refer to [2,18] and bibliography cited therein.
In the present article we adopt the point of view of the continuum approach and we apply such point of view to the numerical study of a single ion channel. This aims at reproducing in mathematical terms the studies that are done in modern Electrophysiology using the patch clamp technique. Our choice is principally motivated by the fact that the scope of the work is to investigate the interplay among several physical mechanisms that are recognized to affect the biological response of ion channels but are far from being assessed in a quantitative manner starting from a first principle analysis. In this respect, we believe that the use of a continuum-based formulation may be convenient from a computational standpoint and sufficiently accurate to provide information about these phenomena.
In the context of formulations based on the continuum assumption, the most widely adopted approach to ion channel modeling and simulation is represented by compartmental models described by a system of ordinary differential equations (ODEs) based on the solution of Kirchhoff current law (KCL) written at the cellular membrane level [27]. The KCL equation is supplemented by phenomenological expressions characterizing the input-output functional response of each protein channel to changes in ion concentrations and/or electric potential inside and outside the cell [19]. Temperature in these models is usually assumed to be a given parameter.
More sophisticated approaches involve the solution of a system of partial differential equations (PDEs) expressing balance of mass of each single ion species flowing across the channel and the Gauss law for the electric field. The system is supplemented by a transport relation, known as the Nernst-Planck (NP) equation, that describes ion motion under an electrochemical gradient [22,24,36]. Also in this kind of modeling, well-known as the Poisson-Nernst-Planck (PNP) system and as the Drift-Diffusion (DD) system [23,30], temperature is a given parameter. A significant step forward to account for temperature as a dependent variable was taken in [12]. In this reference a hydrodynamic (HD) formulation including convective and thermal energy in the electro-chemical motion of a single cation is proposed and numerically investigated. A remarkable feature of [12] is that model and analysis are inspired and guided by the analogy between a biological ion channel and the channel of a semiconductor device in which electrons and holes, instead of charged ion particles, flow to transport electrical current between device terminals. This similarity between biology and solid-state electronics has been thoroughly addressed in the overview paper [14], in the numerical simulations of [17,22,33] [29]. In this latter reference, the mathematical view of [12] is generalized by the introduction of a hierarchical modeling perspective to represent ion transport in a biological channel in which the interstitial electrolyte fluid is assimilated to the semiconductor device medium where two monovalent species (anion and cation) are flowing under the effect of electric, chemical and thermal forces. The hierarchy proposed in [29] is based on the ideas discussed in the semiconductor modeling reference books [41], [31], [23] and [26], and includes four members: the basic DD formulation, the electro-thermal (ET) formulation, the hydrodynamic (HD) formulation and the thermo-hydrodynamic (THD) formulation. For each member of the hierarchy, a set of conservation laws for mass, momentum and energy is written. In the case of DD, ET and HD models, the energy exchange between particles and medium are neglected while in the case of the THD model a supplementary conservation law is added to account for the dynamical thermo-electrochemical balance among the three interacting subsystems. The hierarchy is extensively investigated in a series of numerical computations performed in a simplified one-dimensional channel geometry. Externally applied data on which a sensitivity analysis is carried out are the values of bulk ion concentrations in the intra and extra-cellular sites and the applied potential drop across the channel. Simulations indicate that: i) channel heating is, in general, relatively small compared to ion heating; and ii) for certain ranges of model parameters, the high-order effects introduced by the THD picture can significantly affect the input-output transfer characteristics of the "biological transistor".
Based on the experience and results described above, in the present article we propose the following mathematical structure for the modeling of ion charge transport in biological channels: (a) in Section 2 we review the classic DD (PNP) model for ion electrodiffusion; (b) in Section 3 we use the channel conformation model proposed in [11] to supply the DD (PNP) formulation with a set of "reduced order" boundary conditions that allow one to account for the experimentally accessible values of applied bias, ion concentrations and bathing temperatures; (c) in Section 4 we review two members of the hierarchy of models for ion charge transport introduced in [29], the Thermo-Hydrodynamic and Electro-Thermal systems, and include electroosmotic effects by adding a linear advective term into the momentum balance equations proportional to the given electrolyte fluid velocity. The resulting modified hierarchy is a family of (approximated) velocity-extended models, as the self-consistent study of electrolyte fluid dynamics through Navier-Stokes equations is neglected unlike in [24,25,36,40]; (d) in Section 5 we illustrate the solution map used to solve iteratively the models in (a), (b) and (c) in steadystate conditions; (e) in Section 6 we study a linear model advective-diffusive and reactive boundary value problem (BVP) that represents each of the subproblems arising in the solution map, focusing on the continuous maximum principle; (f) in Section 7 we address the numerical study of the model BVP in (e) based on the dual-mixed hybridized (DMH) finite element method proposed and analyzed in [4,7,34] for elliptic equations and extended to more general differential operators in [3]; (g) in Section 8 we propose a stabilization of the DMH scheme to deal with the case where advection or reaction dominate over diffusion in the model BVP; (h) in Section 9 we conduct a series of numerical tests to verify the convergence and stability properties of the DMH method; (i) in Section 10 we apply external temperature gradients, besides the usual electrochemical gradients, to the biological transistor to replicate in vitro the biophysical situations occurring in vivo in thermo-TRP channels; (j) in the concluding Section 11 we summarize the main contributions of the present research and we address a list of forthcoming activities for model improvement.

2 The Poisson-Nernst-Planck model for ion electrodiffusion
The Poisson-Nernst-Planck system (PNP) for ion electrodiffusion reads [36]: In the PNP system, (1a) is the continuity equation describing mass conservation for each ion whose concentration is denoted by c i (m −3 ), i = 1, . . . , M, M ≥ 1 being the number of ions flowing in the electrolyte fluid. Each ion flux density f i (m −2 s −1 ) is defined by the Nernst-Planck relation (1b) in which it is possible to recognize a chemical contribution and an electric contribution, in such a way that the model be regarded as an extension of Fick's law of diffusion to the case where the diffusing particles are also moved by electrostatic forces with respect to the fluid. The quantity z i is the valence of the i-th ion, while µ i and D i are the mobility and diffusivity of the chemical species, respectively, related by the Einstein relation (1c) where V th = k B Tsys /q (V) is the thermal potential, our having denoted by k B , Tsys and q, the Boltzmann constant, the absolute temperature of the system and the elementary charge, respectively. We remark that in the PNP modeling approach, Tsys is a parameter representing the constant temperature of the system, without distinction between the ionic species and the electrolyte fluid. This assumption is going to be relaxed in Section 4. The function ρ is the space charge density (Cm −3 ) in the electrolyte and is given by the sum of two contributions, the mobile ion charge q ∑︀ M i=1 z i c i and the fixed charge qP, where P = P(x) is a space-dependent charge density profile whose absolute value is measured in (M). The electric field E (Vm −1 ) due to space charge distribution ρ in the electrolyte is determined by the Poisson equation (1d) which represents Gauss' law in differential form, ε being the dielectric permittivity of the electrolyte fluid medium. For further development, it is useful to introduce the electrical current density J i (Am −2 ), equal to the number of ion charges flowing through a given surface area per unit time and defined as Remark 1. The PNP system (1) has the same format and structure as the Drift-Diffusion (DD) equations for semiconductors (see, e.g., [23]), but it is applied to a different medium (water instead of a semiconductor crystal lattice) and includes, in general, more charge carriers than just holes and electrons, as in the case of semiconductor device theory.

Geometry, biophysical assumptions and boundary conditions
Following the presentation of [11], we illustrate in Fig. 1 the schematic representation of a cross-section of a biological channel, assuming rotational invariance around the channel axis (x axis). Five regions can be distinguished, ordered from left to right:  4. Region 4, d ≤ x < d + : this is the channel antichamber from the extracellular side. 5. Region 5, d + ≤ x ≤ R: this is the bathing solution in the extracellular side.
The above geometrical representation leads naturally to a multi-domain formulation for ion channel simulation. The biophysical treatment of such a problem is fully carried out in [11] whilst its mathematical and numerical treatment is the object of the present article, where the PNP formulation (1) is extended to include thermo-hydrodynamical phenomena to ion electrodiffusion in the channel.

Assumptions
The endpoints of the domain, x = L and x = R, are located sufficiently far from the antichamber and channel regions, in such a way that appropriate equilibrium conditions can be applied. More importantly, the endpoints are assumed to be the physical place where the solution intra and extra-cellular electrochemical conditions are accessible to experimental measurements. Because of this, following [10,11], we assume that at x = L and x = R: (A1) the electric potential φ is a known given quantity, so that: (A2) the ion concentrations c i are known given quantities, so that: where M ≥ 1 is the number of ions flowing in the cellular solution. The boundary values for the ion concentrations satisfy the electroneutrality constraint: where z i is the charge number associated with each ion species c i (z i > 0 for cations, z i < 0 for anions and z i = 0 for neutral species).
We close the characterization of the bathing regions by assuming that: (A3) the ion current densities J i (x) vanish inside the baths Continuing to move from the periphery of the domain towards the channel region, we encounter the antichamber openings, at x = 0 − and x = d + , respectively. Located on the membrane lipid bilayer, a spacedependent fixed charge density P = P(x) is distributed. We make the following assumptions within the antichamber regions: (A4) electroneutrality holds at the channel mouth entrances: (A5) the electric potential is constant in the antichamber regions: (2l) (A6) the ion concentrations are constant in the antichamber regions:

Boundary Conditions at Channel Openings
In this section we use the geometrical multi-domain representation of the problem of Sect. 3 and the assumptions made in Sect. 3.1 to derive the boundary conditions to be supplied to the PNP Model at x = 0 and x = d (channel openings). Using Einstein's relation (1c) in (1g) we can write the PNP current density as where where φ ec i,L and φ ec i,R , i = 1, . . . , M, are constants yet to be determined. From (3c) we get Thus, using (3d) and (3e) and assumptions (A5) and (A6), we find: Using (A2) in the equations above at z = L and z = R, respectively, we obtain the boundary values for the electrochemical potential: It is important to note that relations (3g) and (3h) allow one to express the electrochemical potentials at the contacts as a function of the sole accessible quantities. This is in contrast with [29], page 42, relations From the previous discussion, we see that electroneutrality holds at the external contacts (x = L, x = R) and at the channel mouth entrances (x = 0 − , x = d + ), but not, in general, elsewhere. Therefore, it makes sense to introduce the voltage drops occurring in the bathing solution regions, the so-called built-in potentials: Using the charge neutrality conditions (2i) and (2j), the MB relations (3c) and the definitions (4), the builtin potentials can be determined by solving the two following nonlinear algebraic equations: Remark 3. In the particular case of the KCl solution (z 1 = +1, z 2 = −1) considered in [10], equations (4c) and (4d) can be explicitly solved, to yield:  In order to complete the characterization of the boundary values for the dependent variables of the problem we need the value of the electric potential and of the ion concentrations at the two endpoints of the channel. The potential is determined using (4) and (A5), which yield: To determine the ion concentrations we use (3f), (3g) and (4a) at x = 0 to obtain and (3f), (3h) and (4b) at x = d to obtain Summarizing, the Dirichlet boundary conditions for the PNP system in the reduced channel domain Ω = (0, d) are at x = 0: and, at x = d:

Extensions of the PNP model
In this section we present two extensions of the PNP equation system illustrated in Sections 2 and 3. From now on we restrict our attention to the study of a channel in stationary conditions and to simplify the exposition, in close analogy with semiconductor device physics, we consider a binary mixture of monovalent anions and cations, i.e., such that their chemical valence is equal to ±1 as in the case of the KCl solution or the NaCl solution. We use the symbol p (positive) and n (negative) to refer to the concentrations of cations and anions, respectively. The symbol ν is used to indicate either p or n. When the operators ∓ or ± are used, the upper sign always refers to n-type particles, the lower sign to p-type particles, and the value for ν must be chosen accordingly. The quantities associated with the electrolyte medium are referred to by the symbol e, while x denotes henceforth the spatial coordinate.
In our analysis we follow [29], the series of references [5,23,37] in the context of semiconductor device modeling, and the theory of [24,25,36,40]. The first extension is denoted velocity-extended Thermo-Hydrodynamic model (vTHD) and the second extension is denoted velocity-extended Electro-Thermal model (vET). Both extensions include thermal driving forces and the translational contribution of electrolyte fluid flow, besides the usual electrochemical forces accounted for by the PNP equations. In the vTHD model the temperature of anions, cations and lattice are considered as distinct dependent variables, while in the vET model there is only one temperature, generally varying with the spatial coordinate x, to describe the global system.
Finally, in the present article we assume the electrolyte fluid velocity ve to be a given user-defined constant value. This is a strong simplification which is going to be removed in a future publication (cf. Sect. 11). Moreover, to self-consistently account for electrodynamical effects, we adjoin to both vTHD and vET model the solution of the Poisson equation, written below in conservation form: Notice that, in general, the dielectric permittivity is a space dependent function, that is ε = ε(x). The computer code implementing the discretization scheme of Sect. 7 can deal with this situation, as demonstrated in [29] in the numerical study of the K-channel. However, this is not the case of the simulated channels shown in Section 10, where a constant value of ε is prescribed inside the channel region.

Thermodynamical Equilibrium
Let Tp = Tp(x), Tn = Tn(x) and Te = Te(x) denote the temperatures of cations, anions and electrolyte, respectively, with x ∈ [0, d]. Let also Tsys be a given value representing the constant temperature of the whole environment (ions + fluid). By definition, thermal equilibrium is the condition corresponding to applying no external electrical, mechanical, chemical and thermal forces to the biophysical system. In such a case, system response is represented by the following values of the dependent variables: Eq. (8a) expresses the fact that the ion drift velocity is equal to zero, while Eq. (8b) expresses the fact that the fluid translational velocity is equal to zero. Eq. (8c) expresses the fact that the system settles in an isothermal condition.

Definition 1 (Consistency). A model is consistent if under thermal equilibrium (i.e., if conditions (8) hold) the following conditions are satisfied:
where Jn and Jp are anion and cation current densities, while Sn, Sp and Se are ion and fluid energy fluxes.

Remark 4 (Electrochemical potentials at thermal equilibrium). Using (3a) in (9a) tells us that at thermal equilibrium each ion electrochemical potential is constant in
where the constant value φ ec ν,eq is different for each ion and can be computed as detailed in Sect. 3

Velocity-extended Thermo-Hydrodynamic (vTHD) model
The stationary velocity-extended thermo-hydrodynamic model in the one-dimensional setting consists of the following system of conservation laws (see [41,Chap. 2] and [5,23,37]): In system (11), ν and Tν represent ion density and temperature, Ne and Te are electrolyte number density and temperature, while vν, Jν, Sν and Se represent the drift (or translational) velocity of each ion, the ion current density, the ion energy density flux and the electrolyte energy density flux, respectively. The drift velocities and the current densities are related through the phenomenological law Eq. (11a) expresses conservation of ion density, eq. (11b) expresses conservation of ion momentum density, eqns. (11c) and (11d) express conservation of the energy densities for ions and electrolyte fluid, eqns. (11e) and (11f) are constitutive laws for ion and electrolyte energy density fluxes Sν and Se, respectively.
The quantities κν and κe are the thermoconductivity coefficients of ions and electrolyte, respectively, while τp ν and τw ν are the momentum and energy relaxation times, respectively. Phenomenological expressions for these coefficients are: where µ 0ν is the low-field mobility and v sat is the saturation velocity. For a physical interpretation of the parameter v sat and its influence on system behavior, we refer to [12,29]. The ion mobility is related to the diffusion coefficient Dν by the generalized Einstein relation where mν is the ionic mass. Finally, wν and we are the ion and electrolyte energies, in which we can distinguish a kinetic contribution and a thermal contribution. The quantities w eq p and w eq n are the equilibrium energies of cations and anions.  Proof. Thermal equilibrium implies that conditions (8) are satisfied. In particular, using (8a) into (12a) we get Jp(x) = Jn(x) = 0, i.e., conditions (9a). Using (8b) and (8c) in (11b) and (14b) we get or, equivalently, The above equation shows that for the vTHD model at thermal equilibrium the electrochemical potential (10) is satisfied. Having proved (9a) (and, equivalently, (10)) and using (8c) in (11e) and (11f), we immediately obtain that is the energy flux equilibrium condition (9b). We notice that at thermal equilibrium the ion and fluid energies coincide with their rest energy

Velocity-extended Electro-Thermal (vET) model
This model differs from the vTHD because of the following simplifying assumptions: (H1) cations, anions and electrolyte are in local equilibrium, i.e.
This situation occurs when the system is assumed to have had enough time to evolve and to reach steady state; (H2) the kinetic energy is small compared to the thermal energy, i.e.
Thus we can write and neglect the nonlinear inertial term in the momentum balance equation (11b).
Using (H1), and summing up the three resulting energy balance equations we obtain the following equations of the velocity-extended ET system:  (16b), we see that it now furnishes an explicit expression for the current density as a function of concentration, electric potential and temperature. This allows one to adopt the basic ideas of the Scharfetter-Gummel discretization method [39], which is the most widely used numerical scheme in contemporary semiconductor device simulation because of its remarkable stability and accuracy. Remark 6. Simple manipulations of (16b) show that it can be written in the following equivalent DD form The effective electric drift Eν acting on each ion species can be interpreted as the linear superposition of the electric, fluid and thermal forces that drive ion motion in the electrolyte channel fluid. Relation (16f) can be interpreted as the formal limit of (11b) as the relaxation time τp ν → 0.

Remark 7.
The two previous remarks indicate that the vET model differs quite significantly from the vTHD model because the latter is a PDE system with a markedly hyperbolic character due to the convective terms in the momentum balance equations whereas the former is an incompletely parabolic PDE system. We refer the reader to [23] and bibliography cited therein for a deeper analysis of this important issue.
Proposition 2 (Consistency of the vET model). The vET model is consistent in the sense of Def. 1.

Solution map for the vTHD system
In this section we describe the solution map that is used to iteratively solve the vTHD and vET models. The description of the method is discussed in detail in the case of the vTHD system as a similar approach is used to treat also the vET model. The adopted solution map is an extension of the decoupled algorithm known as Gummel's map, a functional tool widely employed in contemporary Drift-Diffusion simulation of semiconductor devices [23,31,41]. The Gummel solution map is basically a nonlinear block Gauss-Seidel iteration that allows one to subdivide the considered PDE model system into blocks of single equations to be solved separately and in sequence until convergence is achieved.
To describe the method, we follow the idea proposed in [1] (cf. Eq. (25) of this latter reference), and we introduce the generalization of the Maxwell-Boltzmann statistics (3c) to the case of the vTHD system: where, for notational brevity, we have omitted the superscript ec in the electrochemical potentials for anions and cations. The main difference between the definitions (18) and the corresponding (3c) valid in the case of the standard PNP model, is that the ion temperatures Tn and Tp are used instead of the system (constant) temperature Tsys. Assuming local thermal equilibrium at the channel entrance and outlet, the temperature of the ions and of the fluid satisfy the following Dirichlet boundary conditions: We notice that T L and T R are the externally accessible temperatures of the intra and extracellular baths, not necessarily being equal to the same value. This feature of the model is important when an external temperature gradient has to be enforced across the channel to model the heat biosensors described in Sect. 1.
The Gummel algorithm for the iterative solution of the vTHD model starts with an initial guess ]︀ T that satisfies the boundary conditions (6) and (19). Then, for k ≥ 0 until convergence, the algorithm consists of the following steps: 1. solve with the damped Newton method (see [41], Chapt. 7) the nonlinear Poisson equation (NLP) resulting from substituting (18) into (7). This returns the updated potential φ k+1 ; 2. solve the continuity and momentum balance equations (11a)-(11b). This returns the updated concentrations p k+1 , n k+1 ; 3. solve the energy equations (11c)-(11e) . This returns the updated ion temperatures T k+1 p , T k+1 n ; 4. solve the fluid energy equation (11d)-(11f). This returns the updated fluid temperature T k+1 e ; 5. update the electrochemical potentials by inverting (18): 6. check the convergence of the iteration by controlling whether the maximum absolute difference between two consecutive iterations k and k + 1 is less than a prescribed tolerance where toll is a given tolerance, U = {φ, φp , φn , Tp , Tn , Te}, and ‖ ·||∞ is the norm in the space L ∞ (Ω). If condition (20c) is satisfied, the algorithm stops.
A deep analysis of the convergence properties and mathematical foundations of the Gummel decoupling algorithm for the stationary DD model of semiconductors can be found in [23]. Preconditioning strategies for improving the convergence rate of the map are investigated in [28].

The Advection-Diffusion-Reaction Model Problem
As outlined in items (e)-(h) of the introduction, we are providing a detailed description of the algorithm, beginning in this section, and continuing through Sect. 9.
All the linearized equations of the vTHD system (11) that need be solved during the iterative procedure described in Sect. 5 In (21), Ω is the open interval (0, d), u and J denote the primal unknown and the associated flux, while c and g are the reaction and production terms such that (21a) can be regarded as a stationary conservation law for u. In view of the numerical treatment of (21) we assume D, v, c and g to be piecewise smooth functions in Ω, the diffusion coefficient D being a positive bounded function while c and g are nonnegative given functions. Dirichlet boundary conditions are expressed by the nonnegative function u : ∂Ω → R + . Let The following properties express important biophysical features of the solution of the linearized BVP (21). We refer the reader to [35] for details and examples.
Definition 2 (Inverse monotonicity). Let w ∈ C 2 (Ω) ∩ C 0 (Ω). We say that L is inverse-monotone if the inequalities Remark 8. Inverse monotonicity expresses the fact that the dependent variable of the problem, say, a concentration, a temperature or a mass density, cannot take negative values.
Theorem 1 (Comparison principle). Suppose that there exists a function ϕ ∈ C 2 (Ω) ∩ C 0 (Ω) such that and we say that ϕ is a barrier function for u.
Combining (22d) and (22g), we obtain the following result which is a very useful tool in the approximation process of the BVP (21).
and we say that u satisfies a continuous maximum principle (CMP).

Finite Element Approximation
The content of the present methodological section is organized as follows. In Sect. 7.1 we introduce the dualmixed weak formulation of (21) and in Sect. 7.2 its corresponding hybridized finite element approximation, denoted DMH method. Then, in Sections 7.3 and 7.4 we illustrate how to reduce the computational complexity of the dual-mixed hybridized method by the use of static condensation, which leads to solving a linear algebraic system of the same structure as a standard nodal-based finite element scheme.

Continuous formulation
In the dual mixed method, the variables u and J in (21) are treated as independent unknowns, each of which belongs to a different functional space, namely u ∈ V := L 2 (Ω) and J ∈ Q := H(div, Ω) ≡ H 1 (Ω). We proceed by formally multiplying equation (21a) by a test function w ∈ V and equation (21b) by a test function q ∈ Q, and integrate by parts over Ω to obtain the problem: Existence and uniqueness of the solution pair (u, J) of problem (23) can be proved under suitable coercivity assumptions, see [13].

Discrete formulation
We where V h is the vector space of piecewise constant polynomials defined over Ω, Q h is the vector space of piecewise linear polynomials discontinuous over Ω, and Λ h,ρ is the space of functions defined only at the vertices of T h with given boundary values. Notice that dim(Q h ) = 2N el because its functions are not continuous at interelement interfaces. The DMH formulation of problem (23) is: where D h : T h → R + is a modified diffusion coefficient containing a stabilization term to deal with advectivedominated problems.

Remark 9.
The proof of existence and uniqueness of the solution (J h , u h , λ h ) of (25) is not an easy task because of the advective term. In what follows, we assume that such a solution always exists and we refer to [7,34] for the general theoretical tools required in the analysis of (25) and to [3] for an example of such analysis in the case of mixed methods interpreted as nonconforming discretizations.

Remark 10. Conditions (25c) imply the continuity of J h across interdomain boundaries.
Remark 11. The hybrid variable λ h is the Lagrange multiplier associated with the continuity constraint of J h · n and it can be regarded as an approximation of u at the grid nodes:

Static condensation
Consider a generic element K ∈ T h . From (25a), we can compute J h | K by inverting the local flux mass matrix as where Now, using (26b) in (26a), we can express J K in terms of the sole hybrid variable λ K and of g K as The above procedure is well known as static condensation and corresponds to Gaussian elimination of the internal variables J h | K and u h | K in favor of λ h | ∂K at the level of the element K ∈ T h (see [4]).

The algebraic system
Enforcing the continuity of J h at each internal node through (25c) yields the following N el − 1 conditions which can be written in the form of a linear tridiagonal algebraic system for λ, i.e., for i = 1, . . . , M − 1: where the entries of the stiffness matrix M ∈ R M−1×M−1 are defined as: and the entries of the load vector The explicit expressions of these entries will be specified in Sect. 8.1.

The Stabilized DMH Method and the Discrete Maximum Principle
If the solution of a given boundary value problem satisfies a CMP, then a properly designed approximation should behave in the same way. A numerical scheme that does not generate spurious global extrema in the interior of the computational domain is said to satisfy a discrete maximum principle (DMP). We assume that the solution u of the BVP (21) satisfies the a priori estimate (22h) and we characterize the conditions under which the same estimate is satisfied also by the DMH approximation λ h computed by solving the linear algebraic system (27c) associated with the problem.
The following definitions turn out to be useful.

Definition 3 (Inverse monotone matrix). An invertible square matrix A of size n is said to be inverse monotone if
the inequality being understood in the element-wise sense.
A special class of monotone matrices is that introduced below.

Definition 4 (M-matrix). An invertible square matrix A is an M-matrix if:
Given a function η h ∈ Λ h,ρ , for any given function ρ : ∂Ω → R, we denote by ρ * h the piecewise linear continuous interpolate of η h over T h .

Theorem 3 (Sufficient condition for DMP). Assume that matrix M is an M-matrix. Then λ
The following (necessary and sufficient) condition is useful to verify the property of being an M-matrix.

Theorem 4 (Discrete comparison principle). Let A be an invertible matrix with non-positive off diagonal entries (A ij ≤ 0 for i ≠ j). Then, A is an M-matrix if and only if there exists a positive vector e such that Ae ≥ 0 (in the component-wise sense)
, with at least one row index i * such that (Ae) i * > 0.

The stabilized DMH method
The solution of the model BVP (21) may exhibit sharp boundary and/or internal layers according to the relative weight of the advective and reactive coefficients with respect to the diffusion term (see [35]). Correspondingly, this is well known to reflect into numerical instability (spurious unphysical oscillations) if the local mesh size h is not sufficiently small. However, a smaller value of h means a larger size of the linear algebraic Brought to you by | Politecnico di Milano Authenticated Download Date | 7/30/19 2:28 PM system and thus an increased computational effort. To avoid this inconvenience, we introduce in the DMH formulation (25) two specific approaches that ensure numerical stability without necessarily resorting to a small value of h: (S1) diagonal lumping of the local mass flux matrix A K (for dominant reaction); (S2) addition of an artificial diffusion (for dominant convection).
In particular, in the remainder of the section we show that, when (S1) and (S2) are used in conjunction, the stiffness matrix M is an M-matrix irrespective of the value of h. This property guarantees that system (27c) is uniquely solvable and that the piecewise linear interpolate of the solution, λ * h , satisfies the DMP and, in particular, the a-priori estimate (25b). We refer to [35] for a thorough discussion of continuous and discrete maximum principles enjoyed by the solution of singularly perturbed BVPs and by their corresponding finite element and finite difference approximations. For ease of presentation, we assume henceforth that problem coefficients D, v, c and g are constant positive quantities and that the mesh size is uniform and given by h = 1/N el .

(S1): Lumping of the local flux mass matrix
The stabilization approach to deal with the case where the reaction coefficient c dominates over the diffusion coefficient D in (21) consists of replacing the local flux mass matrix A K in (26a) with a diagonal matrix̂︀ A K such that: The approach (26) is referred to as mass lumping and is equivalent to using the trapezoidal quadrature rule to compute the integral Using (26) in (25a), when q h = ψ 1 , we get while, when q h = ψ 2 , we similarly get Equations (26c) and (26c) express each local degree of freedom of the flux J h | K as the sum of an average advective flux and a finite difference approximation of the Fick diffusive flux flowing between the boundary of K and its center. The two relations are the DMH analogues of those of the dual-mixed stabilized method of [38].

(S2): Artificial diffusion
The stabilization approach to deal with the case where the advection coefficient v dominates over the diffusion coefficient D in (21) where Pe loc := |v|h 2D is the local Péclet number and Φ : R + → R + is a suitable stabilization function to be chosen in such a way that: The approach (27a) is referred to as artificial diffusion. Two special choices of Φ are: In the case where (27e) is adopted, the resulting stabilized DMH coincides with the classic SG exponentially fitted difference scheme (cf. [39]).
The accuracy of the two stabilized DMH methods, as Thus, the SG formulation is far superior than the upwind method in the limit of a small grid size. However, for a large value of the Pèclet number the two stabilized methods are practically equivalent in terms of the amount of added artificial diffusion (see also [38] and [29]).

Matrix form of the stabilized DMH method
Using superscripts − and + to refer to quantities on the interval ], respectively, let us enforce continuity of J h at every interior node x i (i = 2, . . . , N − 1) Then, combining equations (26c) and (26d) with condition (28a), we get By inspection on (28b) the entries of the stiffness matrix M and of the load vector f become: (28e)

Experimental Validation of the DMH method
In this section we numerically verify the convergence and stability properties of the DMH method.

Convergence Analysis
We consider the model BVP (21) In the numerical approximation of the problem, the mesh size is uniform and equal to h = d/N el , with N el = [10, 20, 40, 80, 160, 320, 640, 1280, 2560] T . Table 1 illustrates the convergence history of the DMH method by reporting various norms for the errors u − u h , Π 0 u − u h and u − λ * h , where Π 0 u is the L 2 projection of u onto the finite element space V h defined in (24a). The discrete maximum norm ‖·‖ ∞,h associated with the triangulation T h is defined for each continuous function η(x) : [0, L] → R as ‖η(x)‖ ∞,h := max Table 2 reports the error J − J h in the L 2 -norm and in the H div ≡ H 1 -norm.
The analysis of the asymptotic convergence orders that are predicted by Tables 1 and 2 show that the computed numerical solutions verify the following error estimates Brought to you by | Politecnico di Milano Authenticated Download Date | 7/30/19 2:28 PM The above results are in excellent agreement with the theoretical convergence rates predicted in the elliptic case in [4] and [7,34].

Reaction-Dominated and Advective-Dominated Regimes
In this section we demonstrate the efficacy of the stabilization techniques proposed in Sect. 8.1 in the study of two model problems, special instances of the BVP (21). For ease of presentation, we set L = 1 and we subdivide the computational domain into N el = 10 uniform intervals of size h = 1/N el = 0.1. We also assume that the coefficients D, v, c and g are constant, with g = 1 and we take homogeneous boundary conditions, u(0) = u(1) = 0 (u = 0).

Diffusion-reaction BVP: mass-lumping
In this case we have v = 0. A barrier function for u is ϕ(x) = 1/c. If the diffusion coefficient is small compared to the reaction term, e.g D = 10 −3 , c = 1, spurious (unphysical) oscillations arise near the boundaries, see Figure 2a. In order to eliminate such oscillations, we adopt the mass-lumping stabilization procedure. The numerical solution λ h obtained with mass-lumping is illustrated in Figure 2b

Diffusion-advection BVP: artificial diffusion
In this case we have c = 0. A barrier function for u is ϕ(x) = x/v. If the diffusion coefficient is small compared to the advective term, e.g D = 5 · 10 −3 , v = 1, spurious (unphysical) oscillations arise in the neighborhood of x = 1 and propagate throughout the entire domain polluting the overall quality of the computed solution, see Figure 3a. Numerical solutions λ h obtained with UP and SG stabilization functions are illustrated in Figures  3b and 3c, respectively. In the case of the SG stabilization function, the computed λ h is nodally exact.

Simulation of biological channels
In this concluding section we carry out a thorough validation of the vET and vTHD models in the study of two different biological channels: 1) the Gramicidin-A channel considered in [12]; and 2) the bipolar nanofluidic diode considered in [45]. The main focus of the simulations is on the current-voltage (IV) characteristics of the channel and on how the IV curves are affected by the boundary conditions, especially the temperature of the two bulk regions, T(0) = T L and T(d) = T R . Also, we aim to investigate the electroosmosis effect, which is accounted for by the electrolyte fluid velocity ve. In all the reported computations, if not otherwise stated, we set N el = 200 and use the SG stabilization function (27e), with the exception of a specific case discussed in Section 10.1. Moreover, in the case of the vTHD model we set the saturation velocity v sat equal to 10 ms −1 . We also assume that the channel has a constant cross-sectional area equal to A, whose value may vary from channel to channel. The current I (A) is thus computed as For further information on the simulated channels, we suggest to consult [29] and the references cited therein.

Gramicidin-A channel (ballistic diode)
This channel is thoroughly analyzed in the work [12] that is here used as a benchmark for the biophysical and numerical assessment of models and methods proposed in the present article. To allow comparison between the results of our models and those of [12] we set ve = 0. Unless otherwise specified, the subsequent figures refer to computations with the THD model, meaning that the outcome of the ET model would appear indistinguishable. The channel length d is equal to 2.5 nm, φ R is always set equal to 0 V and φ L = Vapp, Vapp being the externally applied voltage.

Electrochemical variables
The permanent charge profile of the channel is illustrated in Fig. 4 (left). Because of the negative fixed charge the channel is highly selective to ion flow and attracts positive Na + ions while Cl − ions are mainly repelled. This behavior is confirmed by the computed ion concentration profiles shown in Fig. 4 (right) in the case where a voltage drop φ L − φ R = Vapp = 0.1 V is applied across the channel. We first observe that the cation concentration computed by our algorithm is in excellent agreement with that of [12]. We also notice the presence of strong internal layers in the Na + distribution which are effectively captured by the stabilized DMH discretization without introducing spurious oscillations. In the mentioned reference [12] the hydrodynamic equations are solved using a rather different numerical strategy based on the essentially non oscillatory finite differences with shock capturing introduced in [42,43]. The role of the SG stabilization in the quality of the computed ion concentrations is well documented by the results shown in Fig. 5. In this graph, we see on the left side the cation and anion densities obtained by running the DMH discretization scheme without stabilization (Φ = 0) on an uniform grid of N el = 19 elements. The maximum local Pèclet number is in this case equal to 68.6288 so that the solution exhibits a markedly oscillatory behavior and fails to satisfy the DMP. Conversely, the adoption of the SG stabilization with the same number of elements produces the ion concentrations plotted on the right side of Fig. 5. The computed solution is strictly positive and the internal layers at the interfaces between the channel and highly doped regions are captured with the resolution allowed by the roughness of the grid size. The permanent charge profile and the consequent distribution of ions allows us to regard the Gramicidin-A channel as the biophysical analogue of an electronic "ballistic" diode of type p + -p-p + . This analogy is useful in the interpretation of simulation results and supports the idea that "channels are transistors alive" (cf. [14]). In this respect, the region of the channel x ∈ [0, 0.5] nm corresponds to the Source terminal which has the role of emitting cations into the channel, while the region of the channel x ∈ [2, 2.5] nm corresponds to the Drain terminal which has the role of collecting the cations that have travelled throughout the channel under the action of the thermo-electro-chemical forces. Electric potential and electric field profiles are shown in Fig. 6 (left) and Fig. 6 (right), respectively. Again, we observe that the electric potential and field computed by our algorithm are in excellent agreement with those of [12]. By inspection of the electric potential distribution, we see that the positive pole of the bio-electrical system is located at x = 0, nm while the negative pole is at x = 2.5 nm. Thus, the cations move from left to right, giving rise to accumulation of (fixed) negative charge at the entrance of the channel and positive (mobile) charge at the outlet of the channel. These accumulations correspond to the two strong peaks visible in the electric field distribution. To investigate the effect of different boundary values for the temperature at both entrance and outlet of the channel we let T(0) = T(d) range between 270 K and 330 K. Fig. 7 shows the IV curves for the THD model (left) and the ET model (right), respectively. The externally applied voltage Vapp ranges from 0 to 0.1 V. Temperature influence is visible in both models: the higher the bath temperature the lower the current. This agrees with the biophysical insight that an increase of bath temperature corresponds to a higher thermal energy for the ions and, consequently, to a higher energy dissipation because of frictional effects between particles and fluid. We notice also that the ET model predicts a slightly higher current than the THD model because this latter model accounts for convective energy dissipation within the fluid.    from the ion channel. In the two reservoirs, the temperature profile is linear, because the cation concentration is almost constant and the electric field is very small. Then, ion temperature increases in the channel region according with the fact that particles are accelerated by the electric field from left to right, approximately for x ≥ 0.8 nm (cf. Fig. 6, right). However, the effect of the electric field at the two junctions is quite different in the two sets of thermal boundary conditions. If T(0) > T(d), electric and thermal fields act in the same direction (from left to right) so that as T(0) increases, the thermal flow at the two channel boundaries correspondingly increases. If T(d) > T(0) electric and thermal forces act in opposite directions and this contributes to diminishing the thermal flow at the two channel boundaries. It is interesting to notice that ion heating (i.e., the maximum value of ion temperature inside the channel) is considerably larger in the case where T(d) > T(0), because in this condition ion acceleration due to the electric field greatly dominates over the thermal field along the channel (moving from left to right) so that the ion total energy increases and temperature gets larger. Then, once the ions are injected across the junction between channel and Drain they immediately thermalize to the local energy of the reservoir distribution and approximately cool down to the temperature enforced at x = d. Temperature profiles of the electrolytic fluid are almost linear for both choices of the applied thermal drop, see Fig. 9. This behaviour is to be ascribed to the high value of thermal conductivity of the electrolyte fluid compared to that of the ion fluid. This makes the fluid behave as a perfect sink so that its heating is only passively driven by an external temperature gradient according to Fourier's law (11f) (ve = 0 in this case).

Ion velocity
From ion velocity profiles we can inspect how an externally applied temperature gradient may affect ion motion in the channel.
In Fig. 10  larger spread than in the case of the velocity computed by the THD model. This is to be ascribed to the larger thermal conductivity of the (thermally) unified system composed by anions, cations and fluid. The IV curves computed by the two models in the two distinct thermal sets of boundary conditions reflect what reported about ion velocities: a smaller ion velocity means a smaller current flowing in the channel. The externally applied voltage Vapp ranges from 0 to 0.1 V. Fig. 12 (left) shows the IV curves for the THD model when T(0) > T(d) = 300 K. Notice that the current assumes negative values when Vapp is close to zero and the temperature gradient due to thermal BCs is strong enough to counterbalance the affect of the applied voltage. Consistently with the ion velocity profiles of Fig. 10 (left), the current flowing in the channel is lower than in the case of homogeneous thermal BCs. As T(0) increases, the I-V profile is shifted down with respect to the case of homogeneous BCs (T(0) = 300 K solid black line). The case when T(d) > T(0) = 300 K is shown in Fig. 12 (right). In this situation the IV relationships are shifted up, meaning that as T(d) increases, a higher current flows in the channel, just as predicted from the ion velocity profiles of Fig. 10

IV curves with velocity extended THD and ET models
In the simulations conducted so far, we set ve = 0, which corresponds to neglecting the electroosmosis effect. We have also conducted simulations with ve varying in the range from 0 to 0.01 m s −1 , and we found that it has no appreciable influence on the current flowing in the channel, so that no results are shown in this case. However, different conclusions are drawn in the channel presented in the following section.

Bipolar nanofluidic diode
The nanofluidic channel investigated in the present section (whose detailed representation can be found in [29,45]) is called Bipolar (BP) nanofluidic diode because of its similarity with a p-n electronic diode: as a matter of fact, the BP channel can operate into two different states, depending on the sign of the applied voltage: (1)

Electrochemical variables
The input-output behavior of the BP channel is related to the surface permanent charge profile (cf. Fig. 1). The channel has a negative surface charge along half of its length, while the remaining part has positive surface charge (of the same magnitude), see Fig. 14 (left). Thus, anions and cations carry equal weight to channel behavior. This is confirmed by the symmetric spatial distribution of ion concentrations, both in the case of forward and reverse bias, see Fig. 14 (middle) and Fig. 14 (right), respectively. Notice that ion concentration is much higher in forward bias than reverse bias. Electric potential and electric field profile are shown in Fig. 15. One can see that, in the closed-state, carrier flow is inhibited by the potential barrier at the middle of the channel. Conversely, in the open-state, the potential drop enhances ion flow. The computed IV curves show the on-off trend just described. In particular, we see that if Vapp is negative the current is almost equal to zero, while it increases in a nonlinear manner if Vapp is positive (cf. Fig. 16). The predicted marked on-off function mode of the BP channel and the shape of the I-V curves are also in very good qualitative agreement with data reported in the study of heat biosensors in [9,44,46]. This fact is an encouraging motivation to a further use and calibration of the computational models and methodologies proposed in the present article on biophysical novel applications.

Temperature
The IV curves, when T(0) = T(d) assumes lower or higher values than the reference value of 300 K, are shown in Fig. 16 for the THD model (left) and ET model (right), respectively. The externally applied voltage Vapp ranges from -1 to 1 V. The IV curves computed by the THD model follow the same trend, as a function of the bath temperature, as in Sect. 10.1 (cf. Fig. 7): the higher the bath temperature the lower the current flowing in the channel because of increased frictional effects, see Fig 16 (left). The IV curves computed by the ET model in this type of channel differ remarkably from those computed by the THD model with respect to boundary temperature: the higher the bath temperature, the higher the current, see Fig 16 (right). We also observe that the spread in the value of the maximum channel current as a function of bath temperature is much wider for the THD model than for the ET model. This trend was exactly the opposite in the simulation of the Gramicidin-A channel of Sect. 10.1, as demonstrated by Fig. 7. We point out that the application to the BP channel of the theoretical prediction for channel current given by the ideal diode model (see, e.g., [32]) would yield where I is the total current flowing in the diode and I 0 is the saturation current. Thus, according to (30), an increase of channel system temperature turns out into a decrease in channel current, in accordance with the results computed by the THD model. Despite this preliminary indication in favor of the predictions of the THD model, we feel that the strongly different response of the two models when applied to different channel configurations certainly warrants further investigation and will be the object of a further step of our research activity. Also the temperature profiles of electrolytic fluid differ between the ET model and the THD model. The profiles from ET model have a non linear profile, while those from THD model are linear, see Fig. 17 (left) and Fig. 17 (right), respectively, when T(0) ranges from 270 to 330 K.

Electroosmosis
In this type of channel the influence of the electrolytic fluid velocity is shown in Fig 18 (left) for the vTHD model and in Fig. 18 (right) for the vET model. The selected range for the electrolyte fluid velocity (ve from 0 to 0.01 m/s) is rather arbitrary since our model does not account, at the moment, for a self-consistent coupling between the vTHD system and the NS equations for the fluid. As the fluid velocity has positive sign, electrolytic particles flow from left to right. Accordingly, this additional translational force should enhance the movement of positive ions, while it should reduce that of negative ions, since the two ions move in opposite directions. From the IV curves shown, one can see that the additional driving force due to ve lowers the total current flowing in the channel in both models. Analogous results hold for negative values of ve (not shown). The temperature of the electrolytic fluid velocity for different values of ve is shown in Fig. 19. Temperature is only partially affected by ve. The increase in temperature predicted by the vET model (left) is about 3 K, which is quite remarkable, especially if compared to the results obtained with the vTHD model (middle), where temperature changes are much less significant, instead. The very different prediction of the vTHD formulation is related to the parameter v sat that strongly affects the value of the relaxation time parameter τw ν (cf. (13b)) in the energy balance equation (11c): the larger v sat , the smaller τw ν , which corresponds to almost instantaneous restoration of equilibrium conditions in the electrolyte. Indeed, taking a smaller value of the saturation velocity, for example v sat = 1 m/s, leads to a result that is much closer to the thermal profile predicted by the vET model, see Fig. 19 Figure 19: BP diode. ve ranges from 0 to 0.01 ms −1 . vET model (left) differs greatly from vTHD at v sat = 10 m/s (middle), while it is quite close to vTHD at v sat = 1 m/s (right).

Conclusions and Research Perspectives
In the present article we have proposed and numerically investigated a hierarchy of mathematical models for the simulation of thermal, fluid and electrochemical phenomena in biological transmembrane channels. The hierarchy is an extension of the classic Poisson-Nernst-Planck model for ion electrodiffusion and is conducted along the same lines of thought that have guided the development of the so-called hydrodynamic transport model in the analysis of semiconductor devices. To discretize the proposed models we have devised in the 1D case a robust finite element dual-mixed hybridized method that ensures flux conservation, self-equilibrium and satisfaction of a positivity principle for ion concentrations and temperatures. The numerical scheme has been thoroughly studied in several benchmark problems that demonstrate its accuracy and stability. An appropriate solution map is used to successively solve the nonlinear system of equations arising from model hierarchy and the resulting computational tool has been successfully calibrated and validated in the simulation of two realistic biological channels. Future developments of this study include: -time dependent simulations to describe the response of the channel to externally applied stimuli; -self-consistent coupling of the hierarchy with the solution of the Navier-Stokes equations for the electrolyte fluid; -deeper investigation of the dependence of model predictions on biophysical parameters, for instance, saturation velocity that seems to play a critical role in determining the self-heating effect in the electrolyte fluid; -extension of the numerical scheme to 2D and 3D channel simulation; -the coupling with the Hodgkin-Huxley model for the gating variable kinetics regulating the probability of opening/closure of the channel [15,20,21,27]; -a multi-physics/multi-scale coupling between molecular, continuum and lumped parameter simulations to devise a hierarchy of models characterized by different levels of biophysical accuracy and computational cost.

Conflict of interest:
Author state no conflict of interest.