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.
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.
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.
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.
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.
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.
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 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 . We observe that the errors are sensitive to the choice of the numerical flux coefficient for even degree approximations, but are not for odd degree approximations. Numerical experiments are provided to verify the theoretical results.
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.
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.
In this paper, we consider the high order method for solving the linear transport equations under diffusive scaling and with random inputs. To tackle the randomness in the problem, the stochastic Galerkin method of the generalized polynomial chaos approach has been employed. Besides, the high order implicit-explicit scheme under the micro-macro decomposition framework and the discontinuous Galerkin method have been employed. We provide several numerical experiments to validate the accuracy and the stochastic asymptotic-preserving property.
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.
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.
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.
We consider a one-dimensional reaction-diffusion equation describing single- and two-species population dynamics in an advective environment, based on the modeling frameworks proposed by Lutscher et al. in 2006. We analyze the effect of rate of loss of individuals at both the upstream and downstream boundaries. In the single-species case, we prove the existence of the critical domain size and provide explicit formulas in terms of model parameters. We further derive qualitative properties of the critical domain size and show that, in some cases, the critical domain size is either strictly decreasing over all diffusion rates, or monotonically increasing after first decreasing to a minimum. We also consider competition between species differing only in their diffusion rates. For two species having large diffusion rates, we give a sufficient condition to determine whether the faster or slower diffuser wins the competition. We also briefly discuss applications of these results to competition in species whose spatial niche is affected by shifting isotherms caused by climate change.
Stem cell regeneration is an essential biological process in the maintenance of tissue homeostasis; dysregulation of stem cell regeneration may result in dynamic diseases that show oscillations in cell numbers. Cell heterogeneity and plasticity are necessary for the dynamic equilibrium of tissue homeostasis; however, how these features may affect the oscillatory dynamics of the stem cell regeneration process remains poorly understood. Here, based on a mathematical model of heterogeneous stem cell regeneration that includes cell heterogeneity and random transition of epigenetic states, we study the conditions to have oscillation solutions through bifurcation analysis and numerical simulations. Our results show various model system dynamics with changes in different parameters associated with kinetic rates, cellular heterogeneity, and plasticity. We show that introducing heterogeneity and plasticity to cells can result in oscillation dynamics, as we have seen in the homogeneous stem cell regeneration system. However, increasing the cell heterogeneity and plasticity intends to maintain tissue homeostasis under certain conditions. The current study is an initiatory investigation of how cell heterogeneity and plasticity may affect stem cell regeneration dynamics, and many questions remain to be further studied both biologically and mathematically.
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.
In this paper, we study systems of conservation laws in one space dimension. We prove that for classical solutions in Sobolev spaces , with s > 3/2, the data-to-solution map is not uniformly continuous. Our results apply to all nonlinear scalar conservation laws and to nonlinear hyperbolic systems of two equations.
This paper reviews the adaptive sparse grid discontinuous Galerkin (aSG-DG) method for computing high dimensional partial differential equations (PDEs) and its software implementation. The C++ software package called AdaM-DG, implementing the aSG-DG method, is available on GitHub at https://github.com/JuntaoHuang/adaptive-multiresolution-DG. The package is capable of treating a large class of high dimensional linear and nonlinear PDEs. We review the essential components of the algorithm and the functionality of the software, including the multiwavelets used, assembling of bilinear operators, fast matrix-vector product for data with hierarchical structures. We further demonstrate the performance of the package by reporting the numerical error and the CPU cost for several benchmark tests, including linear transport equations, wave equations, and Hamilton-Jacobi (HJ) equations.
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.
In this paper, we propose a novel Local Macroscopic Conservative (LoMaC) low rank tensor method with discontinuous Galerkin (DG) discretization for the physical and phase spaces for simulating the Vlasov-Poisson (VP) system. The LoMaC property refers to the exact local conservation of macroscopic mass, momentum, and energy at the discrete level. The recently developed LoMaC low rank tensor algorithm (arXiv: 2207.00518) simultaneously evolves the macroscopic conservation laws of mass, momentum, and energy using the kinetic flux vector splitting; then the LoMaC property is realized by projecting the low rank kinetic solution onto a subspace that shares the same macroscopic observables. This paper is a generalization of our previous work, but with DG discretization to take advantage of its compactness and flexibility in handling boundary conditions and its superior accuracy in the long term. The algorithm is developed in a similar fashion as that for a finite difference scheme, by observing that the DG method can be viewed equivalently in a nodal fashion. With the nodal DG method, assuming a tensorized computational grid, one will be able to (i) derive differentiation matrices for different nodal points based on a DG upwind discretization of transport terms, and (ii) define a weighted inner product space based on the nodal DG grid points. The algorithm can be extended to the high dimensional problems by hierarchical Tucker (HT) decomposition of solution tensors and a corresponding conservative projection algorithm. In a similar spirit, the algorithm can be extended to DG methods on nodal points of an unstructured mesh, or to other types of discretization, e.g., the spectral method in velocity direction. Extensive numerical results are performed to showcase the efficacy of the method.
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.
In this paper, we propose a finite volume Hermite weighted essentially non-oscillatory (HWENO) method based on the dimension by dimension framework to solve hyperbolic conservation laws. It can maintain the high accuracy in the smooth region and obtain the high resolution solution when the discontinuity appears, and it is compact which will be good for giving the numerical boundary conditions. Furthermore, it avoids complicated least square procedure when we implement the genuine two dimensional (2D) finite volume HWENO reconstruction, and it can be regarded as a generalization of the one dimensional (1D) HWENO method. Extensive numerical tests are performed to verify the high resolution and high accuracy of the scheme.
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.
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 . 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 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.
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.
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.