Dynamics of Advantageous Mutant Spread in Spatial Death-Birth and Birth-Death Moran Models

The spread of an advantageous mutation through a population is of fundamental interest in population genetics. While the classical Moran model is formulated for a well-mixed population, it has long been recognized that in real-world applications, the population usually has an explicit spatial structure which can significantly influence the dynamics. In the context of cancer initiation in epithelial tissue, several recent works have analyzed the dynamics of advantageous mutant spread on integer lattices, using the biased voter model from particle systems theory. In this spatial version of the Moran model, individuals first reproduce according to their fitness and then replace a neighboring individual. From a biological standpoint, the opposite dynamics, where individuals first die and are then replaced by a neighboring individual according to its fitness, are equally relevant. Here, we investigate this death-birth analogue of the biased voter model. We construct the process mathematically, derive the associated dual process, establish bounds on the survival probability of a single mutant, and prove that the process has an asymptotic shape. We also briefly discuss alternative birth-death and death-birth dynamics, depending on how the mutant fitness advantage affects the dynamics. We show that birth-death and death-birth formulations of the biased voter model are equivalent when fitness affects the former event of each update of the model, whereas the birth-death model is fundamentally different from the death-birth model when fitness affects the latter event.

A High Order Accurate Bound-Preserving Compact Finite Difference Scheme for Two-Dimensional Incompressible Flow

For solving two-dimensional incompressible flow in the vorticity form by the fourth-order compact finite difference scheme and explicit strong stability preserving temporal discretizations,we show that the simple bound-preserving limiter in Li et al.(SIAM J Numer Anal 56:3308–3345,2018)can enforce the strict bounds of the vorticity,if the velocity field satisfies a discrete divergence free constraint.For reducing oscillations,a modified TVB limiter adapted from Cockburn and Shu(SIAM J Numer Anal 31:607–627,1994)is constructed without affecting the bound-preserving property.This bound-preserving finite difference method can be used for any passive convection equation with a divergence free velocity field.

Arbitrary High-Order Fully-Decoupled Numerical Schemes for Phase-Field Models of Two-Phase Incompressible Flows

Due to the coupling between the hydrodynamic equation and the phase-field equation in two-phase incompressible flows,it is desirable to develop efficient and high-order accurate numerical schemes that can decouple these two equations.One popular and efficient strategy is to add an explicit stabilizing term to the convective velocity in the phase-field equation to decouple them.The resulting schemes are only first-order accurate in time,and it seems extremely difficult to generalize the idea of stabilization to the second-order or higher version.In this paper,we employ the spectral deferred correction method to improve the temporal accuracy,based on the first-order decoupled and energy-stable scheme constructed by the stabilization idea.The novelty lies in how the decoupling and linear implicit properties are maintained to improve the efficiency.Within the framework of the spatially discretized local discontinuous Galerkin method,the resulting numerical schemes are fully decoupled,efficient,and high-order accurate in both time and space.Numerical experiments are performed to validate the high-order accuracy and efficiency of the methods for solving phase-field models of two-phase incompressible flows.

High-Order Decoupled and Bound Preserving Local Discontinuous Galerkin Methods for a Class of Chemotaxis Models

In this paper,we explore bound preserving and high-order accurate local discontinuous Galerkin(LDG)schemes to solve a class of chemotaxis models,including the classical Keller-Segel(KS)model and two other density-dependent problems.We use the convex splitting method,the variant energy quadratization method,and the scalar auxiliary variable method coupled with the LDG method to construct first-order temporal accurate schemes based on the gradient flow structure of the models.These semi-implicit schemes are decoupled,energy stable,and can be extended to high accuracy schemes using the semi-implicit spectral deferred correction method.Many bound preserving DG discretizations are only worked on explicit time integration methods and are difficult to get high-order accuracy.To overcome these difficulties,we use the Lagrange multipliers to enforce the implicit or semi-implicit LDG schemes to satisfy the bound constraints at each time step.This bound preserving limiter results in the Karush-Kuhn-Tucker condition,which can be solved by an efficient active set semi-smooth Newton method.Various numerical experiments illustrate the high-order accuracy and the effect of bound preserving.

A Vertex Network Model of Arabidopsis Leaf Growth

Biology provides many examples of complex systems whose properties allow organisms to develop in a highly reproducible,or robust,manner.One such system is the growth and development of flat leaves in Arabidopsis thaliana.This mechanistically challenging process results from multiple inputs including gene interactions,cellular geometry,growth rates,and coordinated cell divisions.To better understand how this complex genetic and cellular information controls leaf growth,we developed a mathematical model of flat leaf production.This two-dimensional model describes the gene interactions in a vertex network of cells which grow and divide according to physical forces and genetic information.Interestingly,the model predicts the presence of an unknown additional factor required for the formation of biologically realistic gene expression domains and iterative cell division.This two-dimensional model will form the basis for future studies into robustness of adaxial-abaxial patterning.

Stability and Time-Step Constraints of Implicit-Explicit Runge-Kutta Methods for the Linearized Korteweg-de Vries Equation

This paper provides a study on the stability and time-step constraints of solving the linearized Korteweg-de Vries(KdV)equation,using implicit-explicit(IMEX)Runge-Kutta(RK)time integration methods combined with either finite difference(FD)or local discontinuous Galerkin(DG)spatial discretization.We analyze the stability of the fully discrete scheme,on a uniform mesh with periodic boundary conditions,using the Fourier method.For the linearized KdV equation,the IMEX schemes are stable under the standard Courant-Friedrichs-Lewy(CFL)conditionτ≤λh.Here,λis the CFL number,τis the time-step size,and h is the spatial mesh size.We study several IMEX schemes and characterize their CFL number as a function ofθ=d/h^(2)with d being the dispersion coefficient,which leads to several interesting observations.We also investigate the asymptotic behaviors of the CFL number for sufficiently refined meshes and derive the necessary conditions for the asymptotic stability of the IMEX-RK methods.Some numerical experiments are provided in the paper to illustrate the performance of IMEX methods under different time-step constraints.

Sparse-Grid Implementation of Fixed-Point Fast Sweeping WENO Schemes for Eikonal Equations

Fixed-point fast sweeping methods are a class of explicit iterative methods developed in the literature to efficiently solve steady-state solutions of hyperbolic partial differential equations(PDEs).As other types of fast sweeping schemes,fixed-point fast sweeping methods use the Gauss-Seidel iterations and alternating sweeping strategy to cover characteristics of hyperbolic PDEs in a certain direction simultaneously in each sweeping order.The resulting iterative schemes have a fast convergence rate to steady-state solutions.Moreover,an advantage of fixed-point fast sweeping methods over other types of fast sweeping methods is that they are explicit and do not involve the inverse operation of any nonlinear local system.Hence,they are robust and flexible,and have been combined with high-order accurate weighted essentially non-oscillatory(WENO)schemes to solve various hyperbolic PDEs in the literature.For multidimensional nonlinear problems,high-order fixed-point fast sweeping WENO methods still require quite a large amount of computational costs.In this technical note,we apply sparse-grid techniques,an effective approximation tool for multidimensional problems,to fixed-point fast sweeping WENO methods for reducing their computational costs.Here,we focus on fixed-point fast sweeping WENO schemes with third-order accuracy(Zhang et al.2006[41]),for solving Eikonal equations,an important class of static Hamilton-Jacobi(H-J)equations.Numerical experiments on solving multidimensional Eikonal equations and a more general static H-J equation are performed to show that the sparse-grid computations of the fixed-point fast sweeping WENO schemes achieve large savings of CPU times on refined meshes,and at the same time maintain comparable accuracy and resolution with those on corresponding regular single grids.

Growth of RB Population in the Conversion Phase of Chlamydia Life Cycle

Upon infecting a host cell,the reticulate body(RB)form of the Chlamydia bacteria simply proliferates by binary fission for an extended period.Available data show only RB units in the infected cells 20 hours post infection(hpi),spanning nearly half way through the development cycle.With data collected every 4 hpi,conversion to the elementary body(EB)form begins abruptly at a rapid rate sometime around 24 hpi.By modeling proliferation and conversion as simple birth and death processes,it has been shown that the optimal strategy for maximizing the total(mean)EB population at host cell lysis time is a bang-bang control qualitatively replicating the observed conversion activities.However,the simple birth and death model for the RB proliferation and conversion to EB deviates in a significant way from the available data on the evolution of the RB population after the onset of RB-to-EB conversion.By working with a more refined model that takes into account a small size threshold eligibility requirement for conversion noted in the available data,we succeed in removing the deficiency of the previous models on the evolution of the RB population without affecting the optimal bang-bang conversion strategy.

Review of Computational Approaches to Optimization Problems in Inhomogeneous Rods and Plates

In this paper,we review computational approaches to optimization problems of inhomogeneous rods and plates.We consider both the optimization of eigenvalues and the localization of eigenfunctions.These problems are motivated by physical problems including the determination of the extremum of the fundamental vibration frequency and the localization of the vibration displacement.We demonstrate how an iterative rearrangement approach and a gradient descent approach with projection can successfully solve these optimization problems under different boundary conditions with different densities given.

Bound-Preserving Discontinuous Galerkin Methods with Modified Patankar Time Integrations for Chemical Reacting Flows

In this paper,we develop bound-preserving discontinuous Galerkin(DG)methods for chemical reactive flows.There are several difficulties in constructing suitable numerical schemes.First of all,the density and internal energy are positive,and the mass fraction of each species is between 0 and 1.Second,due to the rapid reaction rate,the system may contain stiff sources,and the strong-stability-preserving explicit Runge-Kutta method may result in limited time-step sizes.To obtain physically relevant numerical approximations,we apply the bound-preserving technique to the DG methods.Though traditional positivity-preserving techniques can successfully yield positive density,internal energy,and mass fractions,they may not enforce the upper bound 1 of the mass fractions.To solve this problem,we need to(i)make sure the numerical fluxes in the equations of the mass fractions are consistent with that in the equation of the density;(ii)choose conservative time integrations,such that the summation of the mass fractions is preserved.With the above two conditions,the positive mass fractions have summation 1,and then,they are all between 0 and 1.For time discretization,we apply the modified Runge-Kutta/multi-step Patankar methods,which are explicit for the flux while implicit for the source.Such methods can handle stiff sources with relatively large time steps,preserve the positivity of the target variables,and keep the summation of the mass fractions to be 1.Finally,it is not straightforward to combine the bound-preserving DG methods and the Patankar time integrations.The positivity-preserving technique for DG methods requires positive numerical approximations at the cell interfaces,while Patankar methods can keep the positivity of the pre-selected point values of the target variables.To match the degree of freedom,we use polynomials on rectangular meshes for problems in two space dimensions.To evolve in time,we first read the polynomials at the Gaussian points.Then,suitable slope limiters can be applied to enforce the positivity of the solutions at those points,which can be preserved by the Patankar methods,leading to positive updated numerical cell averages.In addition,we use another slope limiter to get positive solutions used for the bound-preserving technique for the flux.Numerical examples are given to demonstrate the good performance of the proposed schemes.

Mathematical Modeling of Cell Polarity Establishment of Budding Yeast

The budding yeast Saccharomyces cerevisiae is a powerful model system for studying the cell polarity establishment.The cell polarization process is regulated by signaling molecules,which are initially distributed in the cytoplasm and then recruited to a proper location on the cell membrane in response to spatial cues or spontaneously.Polarization of these signaling molecules involves complex regulation,so the mathematical models become a useful tool to investigate the mechanism behind the process.In this review,we discuss how mathematical modeling has shed light on different regulations in the cell polarization.We also propose future applications for the mathematical modeling of cell polarization and morphogenesis.

Stability Analysis of Inverse Lax-Wendroff Procedure for a High order Compact Finite Difference Schemes

This paper considers the finite difference(FD)approximations of diffusion operators and the boundary treatments for different boundary conditions.The proposed schemes have the compact form and could achieve arbitrary even order of accuracy.The main idea is to make use of the lower order compact schemes recursively,so as to obtain the high order compact schemes formally.Moreover,the schemes can be implemented efficiently by solving a series of tridiagonal systems recursively or the fast Fourier transform(FFT).With mathematical induction,the eigenvalues of the proposed differencing operators are shown to be bounded away from zero,which indicates the positive definiteness of the operators.To obtain numerical boundary conditions for the high order schemes,the simplified inverse Lax-Wendroff(SILW)procedure is adopted and the stability analysis is performed by the Godunov-Ryabenkii method and the eigenvalue spectrum visualization method.Various numerical experiments are provided to demonstrate the effectiveness and robustness of our algorithms.

Verification and Validation of High-Resolution Inviscid and Viscous Conical Nozzle Flows

Capturing elaborated flow structures and phenomena is required for well-solved numerical flows.The finite difference methods allow simple discretization of mesh and model equations.However,they need simpler meshes,e.g.,rectangular.The inverse Lax-Wendroff(ILW)procedure can handle complex geometries for rectangular meshes.High-resolution and high-order methods can capture elaborated flow structures and phenomena.They also have strong mathematical and physical backgrounds,such as positivity-preserving,jump conditions,and wave propagation concepts.We perceive an effort toward direct numerical simulation,for instance,regarding weighted essentially non-oscillatory(WENO)schemes.Thus,we propose to solve a challenging engineering application without turbulence models.We aim to verify and validate recent high-resolution and high-order methods.To check the solver accuracy,we solved vortex and Couette flows.Then,we solved inviscid and viscous nozzle flows for a conical profile.We employed the finite difference method,positivity-preserving Lax-Friedrichs splitting,high-resolution viscous terms discretization,fifth-order multi-resolution WENO,ILW,and third-order strong stability preserving Runge-Kutta.We showed the solver is high-order and captured elaborated flow structures and phenomena.One can see oblique shocks in both nozzle flows.In the viscous flow,we also captured a free-shock separation,recirculation,entrainment region,Mach disk,and the diamond-shaped pattern of nozzle flows.

Superconvergence of Direct Discontinuous Galerkin Methods:Eigen-structure Analysis Based on Fourier Approach

This paper investigates superconvergence properties of the direct discontinuous Galerkin(DDG)method with interface corrections and the symmetric DDG method for diffusion equations.We apply the Fourier analysis technique to symbolically compute eigenvalues and eigenvectors of the amplification matrices for both DDG methods with different coefficient settings in the numerical fluxes.Based on the eigen-structure analysis,we carry out error estimates of the DDG solutions,which can be decomposed into three parts:(i)dissipation errors of the physically relevant eigenvalue,which grow linearly with the time and are of order 2k for P^(k)(k=2,3)approximations;(ii)projection error from a special projection of the exact solution,which is decreasing over the time and is related to the eigenvector corresponding to the physically relevant eigenvalue;(iii)dissipative errors of non-physically relevant eigenvalues,which decay exponentially with respect to the spatial mesh sizeΔx.We observe that the errors are sensitive to the choice of the numerical flux coefficient for even degree P^(2)approximations,but are not for odd degree P^(3)approximations.Numerical experiments are provided to verify the theoretical results.

A Population Dynamics Approach to the Distribution of Space Debris in Low Earth Orbit

The presence of the debris in the Earth’s orbit poses a significant risk to human activity in outer space.This debris population continues to grow due to ground launches,the loss of external parts from space ships,and uncontrollable collisions between objects.A computationally feasible continuum model for the growth of the debris population and its spatial distribution is therefore critical.Here we propose a diffusion-collision model for the evolution of the debris density in the low-Earth orbit and its dependence on the ground-launch policy.We parametrize this model and test it against data from publicly available object catalogs to examine timescales for the uncontrolled growth.Finally,we consider sensible launch policies and cleanup strategies and how they reduce the future risk of collisions with active satellites or space ships.

Numerical Investigations on the Resonance Errors of Multiscale Discontinuous Galerkin Methods for One-Dimensional Stationary Schrödinger Equation

In this paper,numerical experiments are carried out to investigate the impact of penalty parameters in the numerical traces on the resonance errors of high-order multiscale discontinuous Galerkin(DG)methods(Dong et al.in J Sci Comput 66:321–345,2016;Dong and Wang in J Comput Appl Math 380:1–11,2020)for a one-dimensional stationary Schrödinger equation.Previous work showed that penalty parameters were required to be positive in error analysis,but the methods with zero penalty parameters worked fine in numerical simulations on coarse meshes.In this work,by performing extensive numerical experiments,we discover that zero penalty parameters lead to resonance errors in the multiscale DG methods,and taking positive penalty parameters can effectively reduce resonance errors and make the matrix in the global linear system have better condition numbers.

Stability Analysis and Performance Evaluation of Additive Mixed-Precision Runge-Kutta Methods

Additive Runge-Kutta methods designed for preserving highly accurate solutions in mixed-precision computation were previously proposed and analyzed.These specially designed methods use reduced precision for the implicit computations and full precision for the explicit computations.In this work,we analyze the stability properties of these methods and their sensitivity to the low-precision rounding errors,and demonstrate their performance in terms of accuracy and efficiency.We develop codes in FORTRAN and Julia to solve nonlinear systems of ODEs and PDEs using the mixed-precision additive Runge-Kutta(MP-ARK)methods.The convergence,accuracy,and runtime of these methods are explored.We show that for a given level of accuracy,suitably chosen MP-ARK methods may provide significant reductions in runtime.

A Provable Positivity-Preserving Local Discontinuous Galerkin Method for the Viscous and Resistive MHD Equations

In this paper,we construct a high-order discontinuous Galerkin(DG)method which can preserve the positivity of the density and the pressure for the viscous and resistive magnetohydrodynamics(VRMHD).To control the divergence error in the magnetic field,both the local divergence-free basis and the Godunov source term would be employed for the multi-dimensional VRMHD.Rigorous theoretical analyses are presented for one-dimensional and multi-dimensional DG schemes,respectively,showing that the scheme can maintain the positivity-preserving(PP)property under some CFL conditions when combined with the strong-stability-preserving time discretization.Then,general frameworks are established to construct the PP limiter for arbitrary order of accuracy DG schemes.Numerical tests demonstrate the effectiveness of the proposed schemes.

Exponential Time Differencing Method for a Reaction-Diffusion System with Free Boundary

For reaction-diffusion equations in irregular domains with moving boundaries,the numerical stability constraints from the reaction and diffusion terms often require very restricted time step sizes,while complex geometries may lead to difficulties in the accuracy when discretizing the high-order derivatives on grid points near the boundary.It is very challenging to design numerical methods that can efficiently and accurately handle both difficulties.Applying an implicit scheme may be able to remove the stability constraints on the time step,however,it usually requires solving a large global system of nonlinear equations for each time step,and the computational cost could be significant.Integration factor(IF)or exponential time differencing(ETD)methods are one of the popular methods for temporal partial differential equations(PDEs)among many other methods.In our paper,we couple ETD methods with an embedded boundary method to solve a system of reaction-diffusion equations with complex geometries.In particular,we rewrite all ETD schemes into a linear combination of specificФ-functions and apply one state-of-the-art algorithm to compute the matrix-vector multiplications,which offers significant computational advantages with adaptive Krylov subspaces.In addition,we extend this method by incorporating the level set method to solve the free boundary problem.The accuracy,stability,and efficiency of the developed method are demonstrated by numerical examples.

A Blade Altering Toolbox for Automating Rotor Design Optimization

The Blade Altering Toolbox(BAT)described in this paper is a tool designed for fast reconstruction of an altered blade geometry for design optimization purposes.The BAT algorithm is capable of twisting a given rotor’s angle of attack and stretching the chord length along the span of the rotor.Several test cases were run using the BAT’s algorithm.The BAT code’s twisting,stretching,and mesh reconstruction capabilities proved to be able to handle reasonably large geometric alterations to a provided input rotor geometry.The test examples showed that the toolbox’s algorithm could handle any stretching of the blade’s chord as long as the blade remained within the original bounds of the unaltered mesh.The algorithm appears to fail when the net twist angle applied the geometry exceeds approximately 30 degrees,however this limitation is dependent on the initial geometry and other input parameters.Overall,the algorithm is a very powerful tool for automating a design optimization procedure.