Krylov subspace approximations for the exponential Euler method: Error estimates and the harmonic Ritz approximant

We study Krylov subspace methods for approximating the matrix-function vector product φ(tA)b where φ(z) = [exp(z) - 1]/z. This product arises in the numerical integration of large stiff systems of differential equations by the Exponential Euler Method, where A is the Jacobian matrix of the system. Recently, this method has found application in the simulation of transport phenomena in porous media within mathematical models of wood drying and groundwater flow. We develop an a posteriori upper bound on the Krylov subspace approximation error and provide a new interpretation of a previously published error estimate. This leads to an alternative Krylov approximation to φ(tA)b, the so-called Harmonic Ritz approximant, which we find does not exhibit oscillatory behaviour of the residual error.

Stochastic models for regulatory networks of the genetic toggle switch

Bistability arises within a wide range of biological systems from the λ phage switch in bacteria to cellular signal transduction pathways in mammalian cells. Changes in regulatory mechanisms may result in genetic switching in a bistable system. Recently, more and more experimental evidence in the form of bimodal population distributions indicates that noise plays a very important role in the switching of bistable systems. Although deterministic models have been used for studying the existence of bistability properties under various system conditions, these models cannot realize cell-to-cell fluctuations in genetic switching. However, there is a lag in the development of stochastic models for studying the impact of noise in bistable systems because of the lack of detailed knowledge of biochemical reactions, kinetic rates, and molecular numbers. In this work, we develop a previously undescribed general technique for developing quantitative stochastic models for large-scale genetic regulatory networks by introducing Poisson random variables into deterministic models described by ordinary differential equations. Two stochastic models have been proposed for the genetic toggle switch interfaced with either the SOS signaling pathway or a quorum-sensing signaling pathway, and we have successfully realized experimental results showing bimodal population distributions. Because the introduced stochastic models are based on widely used ordinary differential equation models, the success of this work suggests that this approach is a very promising one for studying noise in large-scale genetic regulatory networks.

Stochastic regularization method for backward Cauchy problem in Banach spaces

We consider a stochastic regularization method for solving the backward Cauchy problem in Banach spaces. An order of convergence is obtained on sourcewise representative elements.

Accelerated leap methods for simulating discrete stochastic chemical kinetics

Biologists are increasingly conscious of the critical role that noise plays in cellular functions such as genetic regulation, often in connection with fluctuations in small numbers of key regulatory molecules. This has inspired the development of models that capture this fundamentally discrete and stochastic nature of cellular biology - most notably the Gillespie stochastic simulation algorithm (SSA). The SSA simulates a temporally homogeneous, discrete-state, continuous-time Markov process, and of course the corresponding probabilities and numbers of each molecular species must all remain positive. While accurately serving this purpose, the SSA can be computationally inefficient due to very small time stepping so faster approximations such as the Poisson and Binomial τ-leap methods have been suggested. This work places these leap methods in the context of numerical methods for the solution of stochastic differential equations (SDEs) driven by Poisson noise. This allows analogues of Euler-Maruyuma, Milstein and even higher order methods to be developed through the Itô-Taylor expansions as well as similar derivative-free Runge-Kutta approaches. Numerical results demonstrate that these novel methods compare favourably with existing techniques for simulating biochemical reactions by more accurately capturing crucial properties such as the mean and variance than existing methods.

An analytical and experimental study of the vibration response of a clamped ribbed plate

An analytical solution is presented in this paper for the vibration response of a ribbed plate clamped on all its boundary edges by employing a travelling wave solution. A clamped ribbed plate test rig is also assembled in this study for the experimental investigation of the ribbed plate response and to provide verification results to the analytical solution. The dynamic characteristics and mode shapes of the ribbed plate are measured and compared to those obtained from the analytical solution and from finite element analysis (FEA). General good agreements are found between the results. Discrepancies between the computational and experimental results at low and high frequencies are also discussed. Explanations are offered in the study to disclose the mechanism causing the discrepancies. The dependency of the dynamic response of the ribbed plate on the distance between the excitation force and the rib is also investigated experimentally. It confirms the findings disclosed in a previous analytical study [T. R. Lin and J. Pan, A closed form solution for the dynamic response of finite ribbed plates. Journal of the Acoustical Society of America 119 (2006) 917-925] that the vibration response of a clamped ribbed plate due to a point force excitation is controlled by the plate stiffness when the source is more than a quarter plate bending wavelength away from the rib and from the plate boundary. The response is largely affected by the rib stiffness when the source location is less than a quarter bending wavelength away from the rib.

An advanced meshless method for time fractional diffusion equation

Recently, because of the new developments in sustainable engineering and renewable energy, which are usually governed by a series of fractional partial differential equations (FPDEs), the numerical modelling and simulation for fractional calculus are attracting more and more attention from researchers. The current dominant numerical method for modeling FPDE is Finite Difference Method (FDM), which is based on a pre-defined grid leading to inherited issues or shortcomings including difficulty in simulation of problems with the complex problem domain and in using irregularly distributed nodes. Because of its distinguished advantages, the meshless method has good potential in simulation of FPDEs. This paper aims to develop an implicit meshless collocation technique for FPDE. The discrete system of FPDEs is obtained by using the meshless shape functions and the meshless collocation formulation. The stability and convergence of this meshless approach are investigated theoretically and numerically. The numerical examples with regular and irregular nodal distributions are used to validate and investigate accuracy and efficiency of the newly developed meshless formulation. It is concluded that the present meshless formulation is very effective for the modeling and simulation of fractional partial differential equations.

Velocity-jump models with crowding effects

Velocity jump processes are discrete random walk models that have many applications including the study of biological and ecological collective motion. In particular, velocity jump models are often used to represent a type of persistent motion, known as a “run and tumble”, which is exhibited by some isolated bacteria cells. All previous velocity jump processes are non-interacting, which means that crowding effects and agent-to-agent interactions are neglected. By neglecting these agent-to-agent interactions, traditional velocity jump models are only applicable to very dilute systems. Our work is motivated by the fact that many applications in cell biology, such as wound healing, cancer invasion and development, often involve tissues that are densely packed with cells where cell-to-cell contact and crowding effects can be important. To describe these kinds of high cell density problems using a velocity jump process we introduce three different classes of crowding interactions into a one-dimensional model. Simulation data and averaging arguments lead to a suite of continuum descriptions of the interacting velocity jump processes. We show that the resulting systems of hyperbolic partial differential equations predict the mean behavior of the stochastic simulations very well.

Models of collective cell motion for cell populations with different aspect ratio : diffusion, proliferation and travelling waves

Continuum, partial differential equation models are often used to describe the collective motion of cell populations, with various types of motility represented by the choice of diffusion coefficient, and cell proliferation captured by the source terms. Previously, the choice of diffusion coefficient has been largely arbitrary, with the decision to choose a particular linear or nonlinear form generally based on calibration arguments rather than making any physical connection with the underlying individual-level properties of the cell motility mechanism. In this work we provide a new link between individual-level models, which account for important cell properties such as varying cell shape and volume exclusion, and population-level partial differential equation models. We work in an exclusion process framework, considering aligned, elongated cells that may occupy more than one lattice site, in order to represent populations of agents with different sizes. Three different idealizations of the individual-level mechanism are proposed, and these are connected to three different partial differential equations, each with a different diffusion coefficient; one linear, one nonlinear and degenerate and one nonlinear and nondegenerate. We test the ability of these three models to predict the population level response of a cell spreading problem for both proliferative and nonproliferative cases. We also explore the potential of our models to predict long time travelling wave invasion rates and extend our results to two dimensional spreading and invasion. Our results show that each model can accurately predict density data for nonproliferative systems, but that only one does so for proliferative systems. Hence great care must be taken to predict density data for with varying cell shape.

Wound healing angiogenesis : the clinical implications of a simple mathematical model

Nonhealing wounds are a major burden for health care systems worldwide. In addition, a patient who suffers from this type of wound usually has a reduced quality of life. While the wound healing process is undoubtedly complex, in this paper we develop a deterministic mathematical model, formulated as a system of partial differential equations, that focusses on an important aspect of successful healing: oxygen supply to the wound bed by a combination of diffusion from the surrounding unwounded tissue and delivery from newly formed blood vessels. While the model equations can be solved numerically, the emphasis here is on the use of asymptotic methods to establish conditions under which new blood vessel growth can be initiated and wound-bed angiogenesis can progress. These conditions are given in terms of key model parameters including the rate of oxygen supply and its rate of consumption in the wound. We use our model to discuss the clinical use of treatments such as hyperbaric oxygen therapy, wound bed debridement, and revascularisation therapy that have the potential to initiate healing in chronic, stalled wounds.

Modelling the interaction of keratinocytes and fibroblasts during normal and abnormal wound healing processes

The crosstalk between fibroblasts and keratinocytes is a vital component of the wound healing process, and involves the activity of a number of growth factors and cytokines. In this work, we develop a mathematical model of this crosstalk in order to elucidate the effects of these interactions on the regeneration of collagen in a wound that heals by second intention. We consider the role of four components that strongly affect this process: transforming growth factor-beta, platelet-derived growth factor, interleukin-1 and keratinocyte growth factor. The impact of this network of interactions on the degradation of an initial fibrin clot, as well as its subsequent replacement by a matrix that is mainly comprised of collagen, is described through an eight-component system of nonlinear partial differential equations. Numerical results, obtained in a two-dimensional domain, highlight key aspects of this multifarious process such as reepithelialisation. The model is shown to reproduce many of the important features of normal wound healing. In addition, we use the model to simulate the treatment of two pathological cases: chronic hypoxia, which can lead to chronic wounds; and prolonged inflammation, which has been shown to lead to hypertrophic scarring. We find that our model predictions are qualitatively in agreement with previously reported observations, and provide an alternative pathway for gaining insight into this complex biological process.

Computationally efficient methods for solving time-variable-order time-space fractional reaction-diffusion equation

Fractional differential equations are becoming more widely accepted as a powerful tool in modelling anomalous diffusion, which is exhibited by various materials and processes. Recently, researchers have suggested that rather than using constant order fractional operators, some processes are more accurately modelled using fractional orders that vary with time and/or space. In this paper we develop computationally efficient techniques for solving time-variable-order time-space fractional reaction-diffusion equations (tsfrde) using the finite difference scheme. We adopt the Coimbra variable order time fractional operator and variable order fractional Laplacian operator in space where both orders are functions of time. Because the fractional operator is nonlocal, it is challenging to efficiently deal with its long range dependence when using classical numerical techniques to solve such equations. The novelty of our method is that the numerical solution of the time-variable-order tsfrde is written in terms of a matrix function vector product at each time step. This product is approximated efficiently by the Lanczos method, which is a powerful iterative technique for approximating the action of a matrix function by projecting onto a Krylov subspace. Furthermore an adaptive preconditioner is constructed that dramatically reduces the size of the required Krylov subspaces and hence the overall computational cost. Numerical examples, including the variable-order fractional Fisher equation, are presented to demonstrate the accuracy and efficiency of the approach.

Double diffusive magneto-convection fluid flow in a strong cross magnetic field with uniform surface heat and mass flux

In this study, magnetohydrodynamic natural convection boundary layer flow of an electrically conducting and viscous incompressible fluid along a heated vertical flat plate with uniform heat and mass flux in the presence of strong cross magnetic field has been investigated. For smooth integrations the boundary layer equations are transformed in to a convenient dimensionless form by using stream function formulation as well as the free variable formulation. The nonsimilar parabolic partial differential equations are integrated numerically for Pr ≪1 that is appropriate for liquid metals against the local Hartmann parameter ξ . Further, asymptotic solutions are obtained near the leading edge using regular perturbation method for smaller values of ξ . Solutions for values of ξ ≫ 1 are also obtained by employing the matched asymptotic technique. The results obtained for small, large and all ξ regimes are examined in terms of shear stress, τw, rate of heat transfer, qw, and rate of mass transfer, mw, for important physical parameter. Attention has been given to the influence of Schmidt number, Sc, buoyancy ratio parameter, N and local Hartmann parameter, ξ on velocity, temperature and concentration distributions and noted that velocity and temperature of the fluid achieve their asymptotic profiles for Sc ≥ 10:0.

A boundary preserving numerical algorithm for the Wright-Fisher model with mutation

The Wright-Fisher model is an Itô stochastic differential equation that was originally introduced to model genetic drift within finite populations and has recently been used as an approximation to ion channel dynamics within cardiac and neuronal cells. While analytic solutions to this equation remain within the interval [0,1], current numerical methods are unable to preserve such boundaries in the approximation. We present a new numerical method that guarantees approximations to a form of Wright-Fisher model, which includes mutation, remain within [0,1] for all time with probability one. Strong convergence of the method is proved and numerical experiments suggest that this new scheme converges with strong order 1/2. Extending this method to a multidimensional case, numerical tests suggest that the algorithm still converges strongly with order 1/2. Finally, numerical solutions obtained using this new method are compared to those obtained using the Euler-Maruyama method where the Wiener increment is resampled to ensure solutions remain within [0,1].