Numerical integration of the master equation in some models of stochastic epidemiology

The processes by which disease spreads in a population of individuals are inherently stochastic. The master equation has proven to be a useful tool for modeling such processes. Unfortunately, solving the master equation analytically is possible only in limited cases (e.g., when the model is linear), and thus numerical procedures or approximation methods must be employed. Available approximation methods, such as the system size expansion method of van Kampen, may fail to provide reliable solutions, whereas current numerical approaches can induce appreciable computational cost. In this paper, we propose a new numerical technique for solving the master equation. Our method is based on a more informative stochastic process than the population process commonly used in the literature. By exploiting the structure of the master equation governing this process, we develop a novel technique for calculating the exact solution of the master equation--up to a desired precision--in certain models of stochastic epidemiology. We demonstrate the potential of our method by solving the master equation associated with the stochastic SIR epidemic model. MATLAB software that implements the methods discussed in this paper is freely available as Supporting Information S1.     

Benchmark solution for the prediction of temperature distributions during radiofrequency ablation of cardiac tissue

Several studies on radiofrequency (RF) ablation are aimed at accurately predicting tissue temperature distributions by numerical solution of the bioheat equation. This paper describes the development of a solution that can serve as a benchmark for subsequent numerical solutions. The solution was obtained using integral transforms and evaluated using a C program. Temperature profiles were generated at various times and for different convection coefficients. In addition, a numerical model was developed using the same assumptions made in obtaining the benchmark solution. Comparison of surface and axial temperature profiles shows that the two solutions match very closely, cross validating the numerical methods used in evaluating both solutions.     

A Traveling Wave Model for Invasion by Precursor and Differentiated Cells

We develop and investigate a continuum model for invasion of a domain by cells that migrate, proliferate and differentiate. The model is applicable to neural crest cell invasion in the developing enteric nervous system, but is presented in general terms and is of broader applicability. Two cell populations are identified and modeled explicitly; a population of precursor cells that migrate and proliferate, and a population of differentiated cells derived from the precursors which have impaired migration and proliferation. The equation describing the precursor cells is based on Fisher’s equation with the addition of a carrying-capacity limited differentiation term. Two variations of the proliferation term are considered and compared. For most parameter values, the model admits a traveling wave solution for each population, both traveling at the same speed. The traveling wave solutions are investigated using perturbation analysis, phase plane methods, and numerical techniques. Analytical and numerical results suggest the existence of two wavespeed selection regimes. Regions of the parameter space are characterized according to existence, shape, and speed of traveling wave solutions. Our observations may be used in conjunction with experimental results to identify key parameters determining the invasion speed for a particular biological system. Furthermore, our results may assist experimentalists in identifying the resource that is limiting proliferation of precursor cells.     

Accurate modelling of unsteady flows in collapsible tubes

The context of this paper is the development of a general and efficient numerical haemodynamic tool to help clinicians and researchers in understanding of physiological flow phenomena. We propose an accurate one-dimensional Runge–Kutta discontinuous Galerkin (RK-DG) method coupled with lumped parameter models for the boundary conditions. The suggested model has already been successfully applied to haemodynamics in arteries and is now extended for the flow in collapsible tubes such as veins. The main difference with cardiovascular simulations is that the flow may become supercritical and elastic jumps may appear with the numerical consequence that scheme may not remain monotone if no limiting procedure is introduced. We show that our second-order RK-DG method equipped with an approximate Roe's Riemann solver and a slope-limiting procedure allows us to capture elastic jumps accurately. Moreover, this paper demonstrates that the complex physics associated with such flows is more accurately modelled than with traditional methods such as finite difference methods or finite volumes. We present various benchmark problems that show the flexibility and applicability of the numerical method. Our solutions are compared with analytical solutions when they are available and with solutions obtained using other numerical methods. Finally, to illustrate the clinical interest, we study the emptying process in a calf vein squeezed by contracting skeletal muscle in a normal and pathological subject. We compare our results with experimental simulations and discuss the sensitivity to parameters of our model.     

一类具有时滞的微生物连续培养数学模型研究

研究了具有强核函数时滞的微生物连续培养数学模型,利用泛函微分方程理论和数值解法得到系统在一定操作条件下存在Hopf分叉以及分叉值随操作参数变化的规律,并研究了Hopf分叉产生的方向及周期解的稳定性,绘制了周期解的图形和相图.该模型定性地描述了实验中的振荡和过渡现象.最后与弱核函数时滞模型、离散时滞模型进行比较,分析了它们对多态、振荡等动态行为的影响。     

Numerical Solution to Generalized Burgers'-Fisher Equation Using Exp-Function Method Hybridized with Heuristic Computation

In this paper, a new heuristic scheme for the approximate solution of the generalized Burgers''-Fisher equation is proposed. The scheme is based on the hybridization of Exp-function method with nature inspired algorithm. The given nonlinear partial differential equation (NPDE) through substitution is converted into a nonlinear ordinary differential equation (NODE). The travelling wave solution is approximated by the Exp-function method with unknown parameters. The unknown parameters are estimated by transforming the NODE into an equivalent global error minimization problem by using a fitness function. The popular genetic algorithm (GA) is used to solve the minimization problem, and to achieve the unknown parameters. The proposed scheme is successfully implemented to solve the generalized Burgers''-Fisher equation. The comparison of numerical results with the exact solutions, and the solutions obtained using some traditional methods, including adomian decomposition method (ADM), homotopy perturbation method (HPM), and optimal homotopy asymptotic method (OHAM), show that the suggested scheme is fairly accurate and viable for solving such problems.     

Random Search Algorithm for Solving the Nonlinear Fredholm Integral Equations of the Second Kind

In this paper, a randomized numerical approach is used to obtain approximate solutions for a class of nonlinear Fredholm integral equations of the second kind. The proposed approach contains two steps: at first, we define a discretized form of the integral equation by quadrature formula methods and solution of this discretized form converges to the exact solution of the integral equation by considering some conditions on the kernel of the integral equation. And then we convert the problem to an optimal control problem by introducing an artificial control function. Following that, in the next step, solution of the discretized form is approximated by a kind of Monte Carlo (MC) random search algorithm. Finally, some examples are given to show the efficiency of the proposed approach.     

Generation of the nervous impulse and periodic oscillations

Usually the models for the excitation and propagation of the nervous impulse are studied either in the space-clamp situation or on a model axon extended on both sides to infinity. Following Fitzhugh in the present paper the release of an impulse train at the axon hillock is studied within the scope of Fitzhugh's BVP model. The existence and stability of periodic oscillations are studied by direct methods, also the relation to Liénard's equation. The exact correspondence between the BVP model and the socalled Nagumo-equation is established. For typical examples the solutions are computed by numerical methods.     

Random currents through nerve membranes

The linear cable equation with uniform Poisson or white noise input current is employed as a model for the voltage across the membrane of a onedimensional nerve cylinder, which may sometimes represent the dendritic tree of a nerve cell. From the Green's function representation of the solutions, the mean, variance and covariance of the voltage are found. At large times, the voltage becomes asymptotically wide-sense stationary and we find the spectral density functions for various cable lengths and boundary conditions. For large frequencies the voltage exhibits “1/f ^3/2 noise”. Using the Fourier series representation of the voltage we study the moments of the firing times for the diffusion model with numerical techniques, employing a simplified threshold criterion. We also simulate the solution of the stochastic cable equation by two different methods in order to estimate the moments and density of the firing time.     

Analysis of tracer experiments: VIII. Integro-differential equation treatment of partly accessible,partly injectable multicompartment systems

An integro-differential equation treatment of multi-compartment systems is developed which permits formal analysis of the incomplete data which is available from partly accessible, partly injectable systems. New transport functions are defined which can be obtained directly from the experimental data. These functions serve to characterize the communication and topology between different accessible compartments and also the reentrant contributions from inaccessible sites. The method gives solutions consistent with those of the differential equation approach when the system is uniformly contiguous and accessible, more complete solutions than those of the integral equation approach when all measured compartments are injectable, and in addition provides complete or partial solutions for certain otherwise analytically intractable systems. Detailed numerical illustrations of the method are given.     

Thermoregulation through skin under variable atmospheric and physiological conditions

Considering three layers of the skin and subcutaneous region, an attempt has been made to obtain the analytical and numerical solutions for temperature distribution of the bioheat equation. The problem is studied under variable physiological parameters and atmospheric conditions. The role of metabolic heat generation, blood mass flow and perspiration in different thicknesses have been noted. The numerical solutions for different parametric values are shown graphically.     

The Influence of Synaptic Weight Distribution on Neuronal Population Dynamics

The manner in which different distributions of synaptic weights onto cortical neurons shape their spiking activity remains open. To characterize a homogeneous neuronal population, we use the master equation for generalized leaky integrate-and-fire neurons with shot-noise synapses. We develop fast semi-analytic numerical methods to solve this equation for either current or conductance synapses, with and without synaptic depression. We show that its solutions match simulations of equivalent neuronal networks better than those of the Fokker-Planck equation and we compute bounds on the network response to non-instantaneous synapses. We apply these methods to study different synaptic weight distributions in feed-forward networks. We characterize the synaptic amplitude distributions using a set of measures, called tail weight numbers, designed to quantify the preponderance of very strong synapses. Even if synaptic amplitude distributions are equated for both the total current and average synaptic weight, distributions with sparse but strong synapses produce higher responses for small inputs, leading to a larger operating range. Furthermore, despite their small number, such synapses enable the network to respond faster and with more stability in the face of external fluctuations.     

Periodic metabolic systems: Oscillations in multiple-loop negative feedback biochemical control networks

Summary For a general multiple loop feedback inhibition system in which the end product can inhibit any or all of the intermediate reactions it is shown that biologically significant behaviour is always confined to a bounded region of reaction space containing a unique equilibrium. By explicit construction of a Liapunov function for the general n dimensional differential equation it is shown that some values of reaction parameters cause the concentration vector to approach the equilibrium asymptotically for all physically realizable initial conditions. As the parameter values change, periodic solutions can appear within the bounded region. Some information about these periodic solutions can be obtained from the Hopf bifurcation theorem. Alternatively, if specific parameter values are known a numerical method can be used to find periodic solutions and determine their stability by locating a zero of the displacement map. The single loop Goodwin oscillator is analysed in detail. The methods are then used to treat an oscillator with two feedback loops and it is found that oscillations are possible even if both Hill coefficients are equal to one.     

Nonlinear cable equations for axons. I. Computations and experiments with internal current injection

Steady-state potential and current distributions resulting from internal injection of current in the squid giant axon have been measured experimentally and also computed from nonlinear membrane cable equation models by numerical methods, using the Hodgkin-Huxley equations to give the membrane current density. The solutions obtained by this method satisfactorily reproduce experimental measurements of the steady-state distribution of membrane potential. Computations of the input current-voltage characteristic for a nonlinear cable were in excellent agreement with measurements on axons. Our results demonstrate the power of Cole's equation to extract the nonlinear membrane characteristics simply from measurement of the input resistance.     

Estimating volume flow rates through xylem conduits

We discuss the errors in common approximations of the volume flow rate for laminar flow through conduits with noncircular transverse sections. Before calculating flow rates, ideal geometric shapes are chosen to represent the noncircular transverse sections. The Hagen–Poiseuille equation used with hydraulic diameter underestimates the volume flow rate for laminar flow through conduits even with such ideal shapes. Correction factors that have been proposed for the Hagen–Poiseuille equation also lead to underestimates of the volume flow rate for those shapes. The exact solutions are sometimes difficult to attain, but rates calculated using the exact solutions for the ideal shapes may be as much as five times higher than the approximated rates for common transversely elongated shapes. Either the exact solutions or the approximations may be used to calculate volume flow rates through the xylem of plants. Both of these methods actually approximate flow through the original conduits because the shapes used are approximations of the conduits’ transverse sections. We recommend using the exact solutions whenever possible; they should be closer to the real solution than other approximations. We give tables of correction factors for use in the cases where calculating volume flow rate from the approximate solution, the Hagen–Poiseuille equation, is more feasible. Obtaining theoretical volume flow rates that are larger than previously thought highlights the need to clarify the causes of differences between the theoretical rates and the smaller measured volume flow rates in plant xylem.     

Complete Numerical Solution of the Diffusion Equation of Random Genetic Drift

A numerical method is presented to solve the diffusion equation for the random genetic drift that occurs at a single unlinked locus with two alleles. The method was designed to conserve probability, and the resulting numerical solution represents a probability distribution whose total probability is unity. We describe solutions of the diffusion equation whose total probability is unity as complete. Thus the numerical method introduced in this work produces complete solutions, and such solutions have the property that whenever fixation and loss can occur, they are automatically included within the solution. This feature demonstrates that the diffusion approximation can describe not only internal allele frequencies, but also the boundary frequencies zero and one. The numerical approach presented here constitutes a single inclusive framework from which to perform calculations for random genetic drift. It has a straightforward implementation, allowing it to be applied to a wide variety of problems, including those with time-dependent parameters, such as changing population sizes. As tests and illustrations of the numerical method, it is used to determine: (i) the probability density and time-dependent probability of fixation for a neutral locus in a population of constant size; (ii) the probability of fixation in the presence of selection; and (iii) the probability of fixation in the presence of selection and demographic change, the latter in the form of a changing population size.     

Numerical integration of population models satisfying conservation laws: NSFD methods

Population models arising in ecology, epidemiology and mathematical biology may involve a conservation law, i.e. the total population is constant. In addition to these cases, other situations may occur for which the total population, asymptotically in time, approach a constant value. Since it is rarely the situation that the equations of motion can be analytically solved to obtain exact solutions, it follows that numerical techniques are needed to provide solutions. However, numerical procedures are only valid if they can reproduce fundamental properties of the differential equations modeling the phenomena of interest. We show that for population models, involving a dynamical conservation law the use of nonstandard finite difference (NSFD) methods allows the construction of discretization schemes such that they are dynamically consistent (DC) with the original differential equations. The paper will briefly discuss the NSFD methodology, the concept of DC, and illustrate their application to specific problems for population models.     

On periodic solutions of a delay integral equation modelling epidemics

Summary A delay-integral equation, proposed by Cooke and Kaplan in [1] as a model of epidemics, is studied. The focus of this work is on the qualitative behavior of solutions as a certain parameter is allowed to vary. It is shown that if a certain threshold is not exceeded then solutions tend to zero exponentially while if this threshold is exceeded, periodic solutions exist. Many features of the numerical studies in [1] are explained.     

Determination of molecular parameters by fitting sedimentation data to finite-element solutions of the Lamm equation.

A method for fitting experimental sedimentation velocity data to finite-element solutions of various models based on the Lamm equation is presented. The method provides initial parameter estimates and guides the user in choosing an appropriate model for the analysis by preprocessing the data with the G(s) method by van Holde and Weischet. For a mixture of multiple solutes in a sample, the method returns the concentrations, the sedimentation (s) and diffusion coefficients (D), and thus the molecular weights (MW) for all solutes, provided the partial specific volumes (v) are known. For nonideal samples displaying concentration-dependent solution behavior, concentration dependency parameters for s(sigma) and D(delta) can be determined. The finite-element solution of the Lamm equation used for this study provides a numerical solution to the differential equation, and does not require empirically adjusted correction terms or any assumptions such as infinitely long cells. Consequently, experimental data from samples that neither clear the meniscus nor exhibit clearly defined plateau absorbances, as well as data from approach-to-equilibrium experiments, can be analyzed with this method with enhanced accuracy when compared to other available methods. The nonlinear least-squares fitting process was accomplished by the use of an adapted version of the "Doesn't Use Derivatives" nonlinear least-squares fitting routine. The effectiveness of the approach is illustrated with experimental data obtained from protein and DNA samples. Where applicable, results are compared to methods utilizing analytical solutions of approximated Lamm equations.