The long-time behavior of the solution to a non-linear diffusion problem in population genetics and combustion

The Fisher (1937) or Kolmogoroff-Petrovsky-Piscounoff (1937) equation exemplifies wave-like phenomena occurring in population genetics and combustion. In an earlier paper, we proposed an extension of this equation and obtained closed form traveling wave, stationary, and “symmetric” solutions. Employing the transformation properties of the extended equation, two integral invariants for the problem are obtained and two Lyapunov functionals, which characterize the evolution of the profile to a uniformly propagating traveling wave, are constructed. A generalization of this modified Fisher equation is proposed and we obtain its integral invariants, traveling wave solutions and wave speeds, as well as the Lyapunov functionals which describe its asymptotic evolution.     

An introduction to the mathematical structure of the Wright–Fisher model of population genetics

In this paper, we develop the mathematical structure of the Wright–Fisher model for evolution of the relative frequencies of two alleles at a diploid locus under random genetic drift in a population of fixed size in its simplest form, that is, without mutation or selection. We establish a new concept of a global solution for the diffusion approximation (Fokker–Planck equation), prove its existence and uniqueness and then show how one can easily derive all the essential properties of this random genetic drift process from our solution. Thus, our solution turns out to be superior to the local solution constructed by Kimura.     

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.     

Can Interactions between Timing of Vaccine-Altered Influenza Pandemic Waves and Seasonality in Influenza Complications Lead to More Severe Outcomes?

Vaccination can delay the peak of a pandemic influenza wave by reducing the number of individuals initially susceptible to influenza infection. Emerging evidence indicates that susceptibility to severe secondary bacterial infections following a primary influenza infection may vary seasonally, with peak susceptibility occurring in winter. Taken together, these two observations suggest that vaccinating to prevent a fall pandemic wave might delay it long enough to inadvertently increase influenza infections in winter, when primary influenza infection is more likely to cause severe outcomes. This could potentially cause a net increase in severe outcomes. Most pandemic models implicitly assume that the probability of severe outcomes does not vary seasonally and hence cannot capture this effect. Here we show that the probability of intensive care unit (ICU) admission per influenza infection in the 2009 H1N1 pandemic followed a seasonal pattern. We combine this with an influenza transmission model to investigate conditions under which a vaccination program could inadvertently shift influenza susceptibility to months where the risk of ICU admission due to influenza is higher. We find that vaccination in advance of a fall pandemic wave can actually increase the number of ICU admissions in situations where antigenic drift is sufficiently rapid or where importation of a cross-reactive strain is possible. Moreover, this effect is stronger for vaccination programs that prevent more primary influenza infections. Sensitivity analysis indicates several mechanisms that may cause this effect. We also find that the predicted number of ICU admissions changes dramatically depending on whether the probability of ICU admission varies seasonally, or whether it is held constant. These results suggest that pandemic planning should explore the potential interactions between seasonally varying susceptibility to severe influenza outcomes and the timing of vaccine-altered pandemic influenza waves.     

Exact moment calculations for genetic models with migration,mutation, and drift

Using properties of moment stationarity we develop exact expressions for the mean and covariance of allele frequencies at a single locus for a set of populations subject to drift, mutation, and migration. Some general results can be obtained even for arbitrary mutation and migration matrices, for example: (1) Under quite general conditions, the mean vector depends only on mutation rates, not on migration rates or the number of populations. (2) Allele frequencies covary among all pairs of populations connected by migration. As a result, the drift, mutation, migration process is not ergodic when any finite number of populations is exchanging genes. In addition, we provide closed-form expressions for the mean and covariance of allele frequencies in Wright's finite-island model of migration under several simple models of mutation, and we show that the correlation in allele frequencies among populations can be very large for realistic rates of mutation unless an enormous number of populations are exchanging genes. As a result, the traditional diffusion approximation provides a poor approximation of the stationary distribution of allele frequencies among populations. Finally, we discuss some implications of our results for measures of population structure based on Wright's F-statistics.     

Diallelic multilocus model of neutral genes

The evolution of a completely linked diallelic multilocus system of neutral genes in a finite population is studied. A diffusion model incorporating random genetic drift and mutation is used. We neglect the recombination. To begin with, the spectral analysis of the Kolmogorov backward equation for this model is investigated. We apply this to two extreme situations when the number of sites approaches to infinity. One is a DeMoivre-Laplace type approximation and the other is a Poisson type approximation. The former is applied to the study of the simultaneous distribution and evolution of a large number of neutral genes. It is applicable to the distribution of a polygenic character controlled by clustered loci on a chromosome, and we show that it differs from the normal distribution on account of random genetic drift and linkage disequilibrium. The latter is applied to the distribution of the number of segregating sites in DNA nucleotide sequences, and the rate of evolution is obtained.     

Analytic computation of the expectation of the linkage disequilibrium coefficient r2

The squared correlation coefficient r(2) (sometimes denoted Delta(2)) is a measure of linkage disequilibrium that is widely used, but computing its expectation E[r(2)] in the population has remained an intriguing open problem. The expectation E[r(2)] is often approximated by the standard linkage deviation sigma(d)(2), which is a ratio of two expectations amenable to analytic computation. In this paper, a method of computing the population-wide E[r(2)] is introduced for a model with recurrent mutation, genetic drift and recombination. The approach is algebraic and is based on the diffusion process approximation. In the limit as the population-scaled recombination rate rho approaches infinity, it is shown rigorously that the asymptotic behavior of E[r(2)] is given by 1/rho+O(rho(-2)), which, incidentally, is the same as that of sigma(d)(2). A computer software that computes E[r(2)] numerically is available upon request.     

Cumulant dynamics of a population under multiplicative selection, mutation, and drift

We revisit the classical population genetics model of a population evolving under multiplicative selection, mutation, and drift. The number of beneficial alleles in a multilocus system can be considered a trait under exponential selection. Equations of motion are derived for the cumulants of the trait distribution in the diffusion limit and under the assumption of linkage equilibrium. Because of the additive nature of cumulants, this reduces to the problem of determining equations of motion for the expected allele distribution cumulants at each locus. The cumulant equations form an infinite dimensional linear system and in an authored appendix Adam Prügel-Bennett provides a closed form expression for these equations. We derive approximate solutions which are shown to describe the dynamics well for a broad range of parameters. In particular, we introduce two approximate analytical solutions: (1) Perturbation theory is used to solve the dynamics for weak selection and arbitrary mutation rate. The resulting expansion for the system's eigenvalues reduces to the known diffusion theory results for the limiting cases with either mutation or selection absent. (2) For low mutation rates we observe a separation of time-scales between the slowest mode and the rest which allows us to develop an approximate analytical solution for the dominant slow mode. The solution is consistent with the perturbation theory result and provides a good approximation for much stronger selection intensities.     

Extinction time and age of an allele in a large finite population

For a single locus with two alleles we study the expected extinction and fixation times of the alleles under the influence of selection and genetic drift. Using a diffusion model we derive asymptotic approximations for these expected exit times for large populations. We consider the case where the fitness of the heterozygote is in between the fitnesses of the homozygotes (incomplete dominance). The asymptotic analysis not only yields new quantitative results but also reveals interesting features that remain hidden in the exact solution. Some of the outcomes are extensions of results known in the literature. The asymptotic approximations also apply to the expected first arrival time of an allele at a specified frequency and to the expected age of an allele.     

A deterministic model for influenza infection with multiple strains and antigenic drift

We describe a multiple strain Susceptible Infected Recovered deterministic model for the spread of an influenza subtype within a population. The model incorporates appearance of new strains due to antigenic drift, and partial immunity to reinfection with related circulating strains. It also includes optional seasonal forcing of the transmission rate of the virus, which allows for comparison between temperate zones and the tropics. Our model is capable of reproducing observed qualitative patterns such as the overall annual outbreaks in the temperate region, a reduced magnitude and an increased frequency of outbreaks in the tropics, and the herald wave phenomenon. Our approach to modelling antigenic drift is novel and further modifications of this model may help improve the understanding of complex influenza dynamics.     

A note on neuronal firing and input variability.

The Ornstein-Uhlenbeck process with a constant forcing function has often been used as a model for the subthreshold membrane potential of a neuron. The mean, variance and coefficient of variation of the first passage time to a constant threshold are examined for this model in the limit of small synaptic noise and low thresholds. A comparison is made between the asymptotic results of Wan & Tuckwell, who used perturbation analysis, and several computationally simpler approximation methods. A generalization of Stein's method gives an overestimate of the mean interval while an approximation by a Wiener process with linear drift gives an underestimate of the mean interval. These bounds are simple to calculate and can be used as a prelude to a more detailed perturbation analysis.     

A cross-diffusion model of forest boundary dynamics

A simple mathematical model of mono-species forest with two age classes which takes into account seed production and dispersal is presented in the paper. This reaction — diffusion type model is then reduced by means of an asymptotic procedure to a lower dimensional reaction — cross-diffusion model. The existence of standing and travelling wave front solutions corresponding to the forest boundary is shown for the later model. On the basis of the analysis, possible changes in forest boundary dynamics caused by antropogenic impacts are discussed.     

Travelling wave solutions of diffusive Lotka-Volterra equations

We establish the existence of travelling wave solutions for two reaction diffusion systems based on the Lotka-Volterra model for predator and prey interactions. For simplicity, we consider only 1 space dimension. The waves are of transition front type, analogous to the travelling wave solutions discussed by Fisher and Kolmogorov et al. for a scalar reaction diffusion equation. The waves discussed here are not necessarily monotone. For any speed c there is a travelling wave solution of transition front type. For one of the systems discussed here, there is a distinguished speed c^* dividing the waves into two types, waves of speed c < c^* being one type, waves of speed c ? c^* being of the other type. We present numerical evidence that for this system the wave of speed c^* is stable, and that c^* is an asymptotic speed of propagation in some sense. For the other system, waves of all speeds are in some sense stable. The proof of existence uses a shooting argument and a Lyapunov function. We also discuss some possible biological implications of the existence of these waves.     

The dynamics of EEG entropy

The scaling properties of human EEG have so far been analyzed predominantly in the framework of detrended fluctuation analysis (DFA). In particular, these studies suggested the existence of power-law correlations in EEG. In DFA, EEG time series are tacitly assumed to be made up of fluctuations, whose scaling behavior reflects neurophysiologically important information and polynomial trends. Even though these trends are physiologically irrelevant, they must be eliminated (detrended) to reliably estimate such measures as Hurst exponent or fractal dimension. Here, we employ the diffusion entropy method to study the scaling behavior of EEG. Unlike DFA, this method does not rely on the assumption of trends superposed on EEG fluctuations. We find that the growth of diffusion entropy of EEG increments of awake subjects with closed eyes is arrested only after approximately 0.5 s. We demonstrate that the salient features of diffusion entropy dynamics of EEG, such as the existence of short-term scaling, asymptotic saturation, and alpha wave modulation, may be faithfully reproduced using a dissipative, first-order, stochastic differential equation—an extension of the Langevin equation. The structure of such a model is utterly different from the “noise+trend” paradigm of DFA. Consequently, we argue that the existence of scaling properties for EEG dynamics is an open question that necessitates further studies.     

A condition for setting off ectopic waves in computational models of excitable cells

The purpose of this paper is to study the stability of steady state solutions of the Monodomain model equipped with Luo-Rudy I kinetics. It is well established that re-entrant arrhythmias can be created in computational models of excitable cells. Such arrhythmias can be initiated by applying an external stimulus that interacts with a partially refractory region, and spawn breaking waves that can eventually generate extremely complex wave patterns commonly referred to as fibrillation. An ectopic wave is one possible stimulus that may initiate fibrillation. Physiologically, it is well known that ectopic waves exist, but the mechanism for initiating ectopic waves in a large collection of cells is poorly understood. In the present paper we consider computational models of collections of excitable cells in one and two spatial dimensions. The cells are modeled by Luo-Rudy I kinetics, and we assume that the spatial dynamics is governed by the Monodomain model. The mathematical analysis is carried out for a reduced model that is known to provide good approximations of the initial phase of solutions of the Luo-Rudy I model. A further simplification is also introduced to motivate and explain the results for the more complicated models. In the analysis the cells are divided into two regions; one region (N) consists of normal cells as model by the standard Luo-Rudy I model, and another region (A) where the cells are automatic in the sense that they would act as pacemaker cells if they where isolated from their surroundings. We let delta denote the spatial diffusion and a denote a characteristic length of the automatic region. It has previously been shown that reducing diffusion or increasing the automatic region enhances ectopic activity. Here we derive a condition for the transition from stable resting state to ectopic wave spread. Under suitable assumptions on the model we provide mathematical and computational arguments indicating that there is a constant eta such that a steady state solution of this system is stable whenever delta approximately > etaa(2), and unstable whenever delta approximately < etaa(2).     

Natural Selection for within-Generation Variance in Offspring Number II. Discrete Haploid Models

In the classical model of genetic drift in population genetics theory, use is made of a hypothetical "infinite-gametic pool". If, instead, the gametic pool is determined by the random number of offspring per individual, a new form of natural selection acting on the variance in offspring number occurs. A diffusion model of this selection process is derived and some of its properties are explored. It is shown that, independent of the sampling scheme used, the diffusion equation has the drift coefficient M(p) = p(1-p) (mul--mu2 + sigma2e2--sigma2el) and the diffusion coefficient v(p) equals p(1-p) [psigma2e2 + (l--p)sigma2el]. It is also pointed out that the Direct Product Branching process model of genetic drift introduces a non-biological interaction between individuals and is thus inappropriate for modeling natural selection.     

Diffusion approximation for a multi-input model neuron

A generalization of an earlier paper (Capocelli and Ricciardi, 1971), dealing with a diffusion approximation for a neuron subject to one excitatory and one inhibitory Poisson input, is provided by not imposing any restrictions on number and magnitude if synaptic inputs. An equation for the neuron's transition p.d.f. is derived, use of which is made to determine the moments of the membrane potential. It is finally shown that a diffusion approximation is possible and that the resulting diffusion process is characterized by constant infinitesimal variance and linear drift.