Solutions to axon equations

The solutions to a general class of axon partial differential equations proposed by FitzHugh which includes the Hodgkin-Huxley equations are studied. It is shown that solutions to the partial differential equations are exactly the solutions to a related set of integral equations. An iterative procedure for constructing the solutions based on standard methods for ordinary differential equations is given and each set of initial values is shown to lead to a unique solution. Continuous dependence of the solutions on the initial values is established and solutions with initial values in a restricted (physiological) range are shown to remain in that range for all time. The iterative procedure is not suggested as the basis for numerical integration.     

Mathematical description of a bursting pacemaker neuron by a modification of the Hodgkin-Huxley equations.

Modifications based on experimental results reported in the literature are made to the Hodgkin-Huxley equations to describe the electrophysiological behavior of the Aplysia abdominal ganglion R15 cell. The system is then further modified to describe the effects with the application of the drug tetrodotoxin (TTX) to the cells' bathing medium. Methods of the qualitative theory of differential equations are used to determine the conditions necessary for such a system of equations to have an oscillatory solution. A model satisfying these conditions is shown to preduct many experimental observations of R15 cell behavior. Numerical solutions are obtained for differential equations satisfying the conditions of the model. These solutions are shown to have a form similar to that of the bursting which is characteristic of this cell, and to preduct many results of experiments conducted on this cell. The physiological implications of the model are discussed.     

Modification of the Kedem-Katchalsky equations

The presented modification of the transport equations of Kedem-Katchalsky resulted in the introduction of (omega s/omega) and omega/(omega-Lp sigma[(1-sigma)C1-(1-sigma s)C2]) factors into the Kedem-Katchalsky equations. The above factors determine the influence of boundary layers on transport across the membrane. The modified Kedem-Katchalsky equations were verified for synthetic membranes and it was shown that the value of the (omega s/omega) factor depended on the type of membrane and the membrane configuration system. This modification facilitated a wider range of application of the Kedem-Katchalsky equations to systems in which the solutions were stirred or unstirred.     

Steady solutions to neoclassical transport equations and the absence of bifurcated states

The question of whether two-valued solutions can exist for an ambipolar electric field in stellarators and rippled tokamaks is considered. Steady solutions to transport equations in the limit of infrequent collisions are obtained in the purely neoclassical transport theory (that is, without allowance for possible anomalous losses). It is shown that, given the particle and heat sources, these equations have only one steady continuous solution, i.e., the steady states are nonbifurcating.     

A gradual update method for simulating the steady-state solution of stiff differential equations in metabolic circuits

Numerical simulation of differential equation systems plays a major role in the understanding of how metabolic network models generate particular cellular functions. On the other hand, the classical and technical problems for stiff differential equations still remain to be solved, while many elegant algorithms have been presented. To relax the stiffness problem, we propose new practical methods: the gradual update of differential-algebraic equations based on gradual application of the steady-state approximation to stiff differential equations, and the gradual update of the initial values in differential-algebraic equations. These empirical methods show a high efficiency for simulating the steady-state solutions for the stiff differential equations that existing solvers alone cannot solve. They are effective in extending the applicability of dynamic simulation to biochemical network models. Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Liquid junction potentials calculated from numerical solutions of the Nernst-Planck and Poisson equations

We present numerical solutions for the one-dimensional Nernst-Planck and Poisson system of equations for steady-state electrodiffusion. Commonly used approximate solutions to these equations invoke assumptions of local electroneutrality (Planck approximation) or constant electric field (Goldman approximation). Calculations were performed to test the ranges over which these approximate theories are valid. For a dilutional junction of a 1:1 electrolyte, separated from adjoining perfectly stirred solutions by sharp boundaries, the Planck approximation is valid for values of kappa dL greater than 10, where 1/kappa d is the Debye length of the more dilute solution. The Goldman approximation is valid for kappa cL less than 0.1 where 1/kappa c is the Debye length of the more concentrated solution. These results suggest that the modeling of electrodiffusive flows in and near membrane ion channels may require numerical solutions of this set of equations rather than the use of either limiting case.     

Ambipolarity condition and possibility of multivalued steady-state solutions to transport equations in stellarators

The possibility of steady-state multivalued solutions to transport equations in stellarators is considered. It is shown that the ambipolarity condition is necessary but not sufficient to find the ambipolar electric field, because the functions entering into it (the plasma density and temperature, as well as their spatial derivatives) depend on the ambipolar field. To do this correctly, it is necessary to solve the full set of time-independent transport equations (including diffusion and heat conduction equations). The possible existence of multivalued solutions to this set of equations is analyzed numerically. It is shown that, under certain conditions that depend on the form and magnitude of particle and heat sources, such solutions can exist. Their form is determined by the initial value of the ambipolar field, the source magnitudes, and the boundary conditions. Discontinuous solutions in which the radial profile of the ambipolar field undergoes jumps are found. In this case, however, the particle and energy fluxes remain continuous, because the discontinuities of the electric field are balanced by the discontinuities of the density and temperature gradients.     

The numerical solution of stiff differential equations

This paper first discusses the conditions in which a set of differential equations should give stable solutions, starting with linear systems assuming that these do not differ greatly in this respect from non-linear systems. Methods of investigating the stability of particular systems are briefly discussed. Most real biochemical systems are known from observation to be stable, but little is known of the regions over which stability persists; moreover, models of biochemical systems may not be stable, because of inaccurate choice of parameter values.The separate problem of stability and accuracy in numerical methods of approximating the solution of systems of non-linear equations is then treated. Stress is laid on the consistently unsatisfactory results given by explicit methods for systems containing "stiff" equations, and implicit multistep methods are particularly recommended for this class of problem, which is likely to include many biochemical model systems. Finally, an iteration procedure likely to give convergence both in multistep methods and in the steady-state approach is recommended, and areas in which improvement in methods is likely to occur are outlined.     

Derivation of general equations describing tracer diffusion in any two-compartment tissue with application to ionic diffusion in cylindrical muscle bundles

A set of generalized diffusion equations have been derived which describe radioactive tracer movement in any tissue that can be modeled as a distributed two-compartment system. These equations have been applied to ionic tracer movement in cylindrical muscle bundles, and the boundary conditions used correspond to experimental conditions during various ionic tracer diffusion experiments on cardiac papillary muscles. Specifically, solutions were obtained for extra- and intracellular tracer washout as well as for the extra- and intracellular steady-state tracer diffusion experiments of Weidmann (1966). These solutions are presented in series form as well as in graphical form and are compared with the corresponding experimental data. A comparison of these solutions with those obtained using simple exponential kinetics is presented, and it is shown that there is a marked discrepancy between these two methods of analysis for bundles of any appreciable diameter.     

Traveling Wave Solutions of a Nerve Conduction Equation

We consider a pair of differential equations whose solutions exhibit the qualitative properties of nerve conduction, yet which are simple enough to be solved exactly and explicitly. The equations are of the FitzHugh-Nagumo type, with a piecewise linear nonlinearity, and they contain two parameters. All the pulse and periodic solutions, and their propagation speeds, are found for these equations, and the stability of the solutions is analyzed. For certain parameter values, there are two different pulse-shaped waves with different propagation speeds. The slower pulse is shown to be unstable and the faster one to be stable, confirming conjectures which have been made before for other nerve conduction equations. Two periodic waves, representing trains of propagated impulses, are also found for each period greater than some minimum which depends on the parameters. The slower train is unstable and the faster one is usually stable, although in some cases both are unstable.     

Accurate computation of the stable solitary wave for the FitzHugh-Nagumo equations

Comparisons are made between three different methods for computing the stable solitary wave solution for the FitzHugh-Nagumo equations which consist of a nonlinear diffusion equation coupled to an ordinary differential equation in time. They model the Hodgkin-Huxley equations which describe the propagation of the nerve impulse down the axon. Two of the methods involve the travelling wave equations. Previous accurate numerical computations of these equations as an initial-value problem using a shooting method lead to inaccurate values for the wave speed; however, nonlinear corrections to the initial values are shown to yield accurate values. A boundary-value method applies asymptotic boundary conditions and uses a spline-collocation code called COLSYS for numerical solution of boundary-value problems which leads to accurate wave profiles and speeds. The third method is to solve an initial-boundary-value problem with an adaptive outgoing wave condition for the partial differential equations where the solitary wave emerges as the stable long time solution. The concept of a wave integral is introduced and they are derived to determine the wave speed used in the adaptive boundary condition and to measure the closeness of the computed solutions to the exact solitary wave solution.This work was supported in part by the Natural Sciences and Engineering Research Council Canada under Grant A4559 and by the John Simon Guggenheim Memorial Foundation     

Isolated periodic solutions of the Hodgkin-Huxley equations

Periodic solutions of the current clamped Hodgkin-Huxley equations (Hodgkin & Huxley, 1952 J. Physiol. 117, 500) that arise by degenerate Hopf bifurcation were studied recently by Labouriau (1985 SIAM J. Math. Anal. 16, 1121, 1987 Degenerate Hopf Bifurcation and Nerve Impulse (Part II), in press). Two parameters, temperature T and sodium conductance gNa were varied from the original values obtained by Hodgkin & Huxley. Labouriau's work proved the existence of small amplitude periodic solution branches that do not connect locally to the stationary solution branch, and had not been previously computed. In this paper we compute these solution branches globally. We find families of isolas of periodic solutions (i.e. branches not connected to the stationary branch). For values of gNa in the range measured by Hodgkin & Huxley, and for physically reasonable temperatures, there are isolas containing orbitally asymptotically stable solutions. The presence of isolas of periodic solutions suggests that in certain current space clamped membrane experiments, action potentials could be observed even though the stationary state is stable for all current stimuli. Once produced, such action potentials will disappear suddenly if the current stimulus is either increased or decreased past certain values. Under some conditions, "jumping" between action potentials of different amplitudes might be observed.     

Exact solutions of Lotka-Volterra equations

The Lotka-Volterra nonlinear differential equations for two competing species P and Q contain six independent parameters. Their general analytic solutions, valid for arbitrary values of the parameters, are at present unknown. However, when two or more of these parameters are interrelated, it is possible to obtain the exact solutions in the P, Q phase plane, and six cases of solvability are given in this paper. The dependence of the solutions on the parameters and the initial conditions can thus be readily investigated.     

Hopf bifurcations in multiple-parameter space of the Hodgkin-Huxley equations I. Global organization of bistable periodic solutions

 The Hodgkin-Huxley equations (HH) are parameterized by a number of parameters and shows a variety of qualitatively different behaviors depending on the parameter values. We explored the dynamics of the HH for a wide range of parameter values in the multiple-parameter space, that is, we examined the global structure of bifurcations of the HH. Results are summarized in various two-parameter bifurcation diagrams with I _ext (externally applied DC current) as the abscissa and one of the other parameters as the ordinate. In each diagram, the parameter plane was divided into several regions according to the qualitative behavior of the equations. In particular, we focused on periodic solutions emerging via Hopf bifurcations and identified parameter regions in which either two stable periodic solutions with different amplitudes and periods and a stable equilibrium point or two stable periodic solutions coexist. Global analysis of the bifurcation structure suggested that generation of these regions is associated with degenerate Hopf bifurcations. Received: 23 April 1999 / Accepted in revised form: 24 September 1999     

A simple vector implementation of the Laplace-transformed cable equations in passive dendritic trees

Transient potentials in dendritic trees can be calculated by approximating the dendrite by a set of connected cylinders. The profiles for the currents and potentials in the whole system can then be obtained by imposing the proper boundary conditions and calculating these profiles along each individual cylinder. An elegant implementation of this method has been described by Holmes (1986), and is based on the Laplace transform of the cable equation. By calculating the currents and potentials only at the ends of the cylinders, the whole system of connected cylinders can be described by a set of n equations, where n denotes the number of internal and external nodes (points of connection and endpoints of the cylinders). The present study shows that the set of equations can be formulated by a simple vector equation which is essentially a generalization of Ohm's law for the whole system. The current and potential n-vectors are coupled by a n × n conductance matrix whose structure immediately reflects the connectivity pattern of the connected cylinders. The vector equation accounts for conductances, associated with driving potentials, which may be local or distributed over the membrane. It is shown that the vector equation can easily be adapted for the calculation of transients over a period in which stepwise changes in system parameters have occurred. In this adaptation it is assumed that the initial conditions for the potential profiles at the start of a new period after a stepwise change can be approximated by steady-state solutions. The vector representation of the Laplace-transformed equations is attractive because of its simplicity and because the structure of the conductance matrix directly corresponds to the connectivity pattern of the dendritic tree. Therefore it will facilitate the automatic generation of the equations once the geometry of the branching structure is known.     

Neural rate equations for bursting dynamics derived from conductance-based equations

A method of obtaining rate equations from conductance-based equations is developed and applied to fast-spiking and bursting neocortical neurons. It involves splitting systems of conductance-based equations into fast and slow subsystems, and averaging the effects of fast terms that drive the slowly varying quantities by showing that their average is closely proportional to the firing rate. The dependence of the firing rate on the injected current is then approximated in the analysis. The resulting behavior of the slow variables is then substituted back into the fast equations, with the further approximation of replacing the fast voltages in these terms by effective values. For bursting neurons the method yields two coupled limit-cycle oscillators: a self-exciting oscillator for the slow variables that commences limit-cycle oscillations at a critical current and modulates a fast spike-generating oscillator, thereby leading to slowly modulated bursts with a group of spikes in each burst. The dynamics of these coupled oscillators are then verified against those of the conductance-based equations. Finally, it is shown how to place the results in a form suitable for use in mean-field equations for neural population dynamics.     

Loss of synchronization in partially coupled Hodgkin-Huxley equations

We study the loss of synchronization of two partially coupled space-clamped Hodgkin-Huxley equations, with symmetric coupling. This models the coupling of two cells through an electrical synapse. For strong enough coupling it is known that all solutions of the equations approach a state where the two cells are perfectly synchronized, having the same behaviour at each moment. We describe the local bifurcations that arise when the coupling strength is reduced, using a mixture of analytical and numerical methods. We find that perfect synchrony is retained for very small positive values of the coupling strength, for almost all initial conditions. Although perfect synchrony is lost for negative values of the coupling constant, the system always retains some degree of synchronization until it becomes totally unstable. This happens in two ways: in many cases for almost all initial conditions the solutions still approach a perfectly synchronized state. Even when this is not true, the attracting solutions are still synchronized, with a half-period phase shift.     

Study of models of single populations: development of equations used in microorganism cultures

This article derives a number of equations which can be used in both continuous and the semicontinuous cultures of microorganism populations in chemostat systems. Using these equations, some phenomena which have been known for many years can be explained reasonably in terms of chemical kinetics, and a number of analytical solutions can be obtained instead of numerical solutions previously published.     

Evolution of discontinuous solutions to transport equations in stellarators

It was shown earlier that, in the range of rare collisions, transport equations for stellarators allow steady discontinuous solutions for the ambipolar electric field and for the plasma density and temperature gradients. Moreover, such solutions are non-single-valued; that is, their explicit form depends on the initial values of the ambipolar electric field. The time-independent transport equations are derived under the conventional quasineutrality condition; i.e., it is assumed that the electron and ion densities, N _e and N _i , are related by the relationship N _e = ZN _i (where Z is the ion charge number). In other words, the plasma charge density is assumed to be much less than the product e _i N _i . Under typical conditions, the corresponding inequality is satisfied by a large margin. However, if the electric field E has discontinuities, then it can be seen from the equation ▿·E = 4πρ that, at the discontinuity points, the charge density becomes infinite and the relationship N _e = ZN _i fails to hold, so it is necessary to replace it with N _e = ZN _i + ρ/e _e . In the transport equations, this latter replacement produces additional terms, proportional to the second radial derivative of the field E. With these additional terms, the steady solutions are modified substantially. First, the ambipolar field and the derivatives of the density and temperatures all become continuous functions of the coordinates, a result that seems to be quite obvious. The second, not-so-obvious result is that the steady solutions become single-valued, i.e., independent of the initial values of the ambipolar electric field. It turns out that, in this case, two regimes are possible, depending on the values of the plasma parameters. In the first regime, the solution is unique and is independent of the initial conditions. In the second regime, two steady solutions can exist, depending on the initial conditions. One of the solution is similar to that obtained in the first regime, and the other differs from the first one both in the ambipolar field profile and in the dependence of the density and temperatures on the minor plasma radius. It cannot be excluded that different plasma confinement modes revealed in experiments are associated with the existence of such solutions.