Search

Weakly nonlinear amplitude equations are derived for the onset of spatially extended patterns on a general class of $n$-component bulk-surface reaction-diffusion systems in a ball, under the assumption of linear kinetics in the bulk. Linear analysis shows conditions under which various pattern modes can become unstable to either generalised pitchfork or transcritical bifurcations depending on the parity of the spatial wavenumber. Weakly nonlinear analysis is used to derive general expressions for the multi-component amplitude equations of different patterned states. These reduced-order systems are found to agree with prior normal forms for pattern formation bifurcations with $O(3)$ symmetry and provide information on the stability of bifurcating patterns of different symmetry types. The analysis is complemented with numerical results using a dedicated finite-element method. The theory is illustrated in two examples; a bulk-surface version of the Brusselator, and a four-component cell-polarity model.

The bulk-surface wave pinning model is a reaction-diffusion system for studying cell polarisation. It is constituted by a surface reaction-diffusion equation, coupled to a bulk diffusion equation with a non-linear boundary condition. Cell polarisation arises as the surface component develops specific patterns. Since proteins diffuse much faster in the cell interior than on the membrane, in the literature, the bulk component is often assumed to be spatially homogeneous. {\re Therefore, the model can be reduced to a single surface equation}. However, {\re in real applications} a spatially non-uniform bulk component might be an important player to take into account. In this paper, we study, through numerical computations, the role of the bulk component and, more specifically, how different bulk diffusion rates might affect the polarisation response. We find that the bulk component is indeed a key factor in {\re determining the surface polarisation response}. Moreover, for certain geometries, it is the spatial heterogeneity of the bulk component that triggers the polarisation response, which might not be possible in a reduced model. Understanding how polarisation depends on bulk diffusivity might be crucial when studying models of migrating cells, which are naturally subject to domain deformation., Comment: 28 pages (last 2 for SM) with 11 figures (2 in the SM)

In this work we present a novel bulk-surface virtual element method (BSVEM) for the numerical approximation of elliptic bulk-surface partial differential equations (BSPDEs) in three space dimensions. The BSVEM is based on the discretisation of the bulk domain into polyhedral elements with arbitrarily many faces. The polyhedral approximation of the bulk induces a polygonal approximation of the surface. Firstly, we present a geometric error analysis of bulk-surface polyhedral meshes independent of the numerical method. Then, we show that BSVEM has optimal second-order convergence in space, provided the exact solution is $H^{2+3/4}$ in the bulk and $H^2$ on the surface, where the additional $\frac{3}{4}$ is due to the combined effect of surface curvature and polyhedral elements close to the boundary. We show that general polyhedra can be exploited to reduce the computational time of the matrix assembly. To demonstrate optimal convergence results, a numerical example is presented on the unit sphere., Comment: 25 pages, 4 figures, 1 table. This replacement improves figures, updates references, and avoids redundancies. arXiv admin note: substantial text overlap with arXiv:2002.11748

Compartmental models are popular in the mathematics of epidemiology for their simplicity and wide range of applications. Although they are typically solved as initial value problems for a system of ordinary differential equations, the observed data is typically akin of a boundary value type problem: we observe some of the dependent variables at given times, but we do not know the initial conditions. In this paper, we reformulate the classical Susceptible-Infectious-Recovered system in terms of the number of detected positive infected cases at different times, we then prove the existence and uniqueness of a solution to the derived boundary value problem and then present a numerical algorithm to approximate the solution.

The mathematical interpretation of interventions for the mitigation of epidemics and pandemics in the literature often involves finding the optimal time to initiate an intervention and/or the use of infections to manage impact. Whilst these methods may work in theory, in order to implement they may require information which is likely not available whilst one is in the midst of an epidemic, or they may require impeccable data about infection levels in the community. In practice, testing and cases data is only as good as the policy of implementation and the compliance of the individuals, which means that understanding the levels of infections becomes difficult or complicated from the data that is provided. In this paper, we aim to develop a different approach to the mathematical modelling of interventions, not based on optimality, but based on demand and capacity of local authorities who have to deal with the epidemic on a day to day basis. In particular, we use data-driven modelling to calibrate an Susceptible Exposed Infectious Recovered-Died (SEIR-D) model to infer parameters that depict the dynamics of the epidemic in a region of the UK. We use the calibrated parameters for forecasting scenarios and understand, given a maximum capacity of hospital healthcare services, how the timing of interventions, severity of interventions, and conditions for the releasing of interventions affect the overall epidemic-picture.

In this paper we consider a coupled bulk-surface PDE in two space dimensions. The model consists of a PDE in the bulk that is coupled to another PDE on the surface through general nonlinear boundary conditions. For such a system we propose a novel method, based on coupling a virtual element method [Beir\~ao da Veiga et al., 2013] in the bulk domain to a surface finite element method [Dziuk & Elliott, 2013] on the surface. The proposed method, which we coin the Bulk-Surface Virtual Element Method (BSVEM) includes, as a special case, the bulk-surface finite element method (BSFEM) on triangular meshes [Madzvamuse & Chung, 2016]. The method exhibits second-order convergence in space, provided the exact solution is $H^{2+1/4}$ in the bulk and $H^2$ on the surface, where the additional $\frac{1}{4}$ is required only in the simultaneous presence of surface curvature and non-triangular elements. Two novel techniques introduced in our analysis are (i) an $L^2$-preserving inverse trace operator for the analysis of boundary conditions and (ii) the Sobolev extension as a replacement of the lifting operator [Elliott & Ranner, 2013] for sufficiently smooth exact solutions. The generality of the polygonal mesh can be exploited to optimize the computational time of matrix assembly. The method takes an optimised matrix-vector form that also simplifies the known special case of BSFEM on triangular meshes [Madzvamuse & Chung, 2016]. Three numerical examples illustrate our findings., Comment: 39 pages, 9 figures, 5 tables

The main purpose of the current paper is to contribute towards the comprehension of the dynamics of the shadow system of a singular Gierer-Meinhardt model on an isotropically evolving domain. In the case where the inhibitor's response to the activator's growth is rather weak, then the shadow system of the Gierer-Meinhardt model is reduced to a single though non-local equation whose dynamics is thoroughly investigated throughout the manuscript. The main focus is on the derivation of blow-up results for this non-local equation, which can be interpreted as instability patterns of the shadow system. In particular, a {\it diffusion-driven instability (DDI)}, or {\it Turing instability}, in the neighbourhood of a constant stationary solution, which then is destabilised via diffusion-driven blow-up, is observed. The latter indicates the formation of some unstable patterns, whilst some stability results of global-in-time solutions towards non-constant steady states guarantee the occurrence of some stable patterns. Most of the derived results are confirmed numerically and also compared with the ones in the case of a stationary domain.

This work presents the development, analysis and numerical simulations of a biophysical model for 3D cell deformation and movement, which couples biochemical reactions and biomechanical forces. We propose a mechanobiochemical model which considers the actin filament network as a viscoelastic and contractile gel. The mechanical properties are modelled by a force balancing equation for the displacements, the pressure and concentration forces are driven by actin and myosin dynamics, and these are in turn modelled by a system of reaction-diffusion equations on a moving cell domain. The biophysical model consists of highly non-linear partial differential equations whose analytical solutions are intractable. To obtain approximate solutions to the model system, we employ the moving grid finite element method. The numerical results are supported by linear stability theoretical results close to bifurcation points during the early stages of cell migration. Numerical simulations exhibited show both simple and complex cell deformations in 3-dimensions that include cell expansion, cell protrusion and cell contraction. The computational framework presented here sets a strong foundation that allows to study more complex and experimentally driven reaction-kinetics involving actin, myosin and other molecular species that play an important role in cell movement and deformation.

Several cellular activities, such as directed cell migration, are coordinated by an intricate network of biochemical reactions which lead to a polarised state of the cell, in which cellular symmetry is broken, causing the cell to have a well defined front and back. Recent work on balancing biological complexity with mathematical tractability resulted in the proposal and formulation of a famous minimal model for cell polarisation, known as the wave pinning model. In this study, we present a three-dimensional generalisation of this mathematical framework through the maturing theory of coupled bulk-surface semilinear partial differential equations in which protein compartmentalisation becomes natural. We show how a local perturbation over the surface can trigger propagating reactions, eventually stopped in a stable profile by the interplay with the bulk component. We describe the behavior of the model through asymptotic and local perturbation analysis, in which the role of the geometry is investigated. The bulk-surface finite element method is used to generate numerical simulations over simple and complex geometries, which confirm our analysis, showing pattern formation due to propagation and pinning dynamics. The generality of our mathematical and computational framework allows to study more complex biochemical reactions and biomechanical properties associated with cell polarisation in multi-dimensions.

This work explores the influence of domain size of a non-compact two dimensional annular domain on the evolution of pattern formation that is modelled by an \textit{activator-depleted} reaction-diffusion system. A closed form expression is derived for the spectrum of Laplace operator on the domain satisfying a set of homogeneous conditions of Neumann type both at inner and outer boundaries. The closed form solution is numerically verified using the spectral method on polar coordinates. The bifurcation analysis of \textit{activator-depleted} reaction-diffusion system is conducted on the admissible parameter space under the influence of two bounds on the parameter denoting the thickness of the annular region. The admissibility of Hopf and transcritical bifurcations is proven conditional on the domain size satisfying a lower bound in terms of reaction-diffusion parameters. The admissible parameter space is partitioned under the proposed conditions corresponding to each case, and in turn such conditions are numerically verified by applying a method of polynomials on a quadrilateral mesh. Finally, the full system is numerically simulated on a two dimensional annular region using the standard Galerkin finite element method to verify the influence of the analytically derived conditions on the domain size for all types of admissible bifurcations., Comment: arXiv admin note: substantial text overlap with arXiv:1801.03151

The ability to locally degrade the extracellular matrix (ECM) and interact with the tumour microenvironment is a key process distinguishing cancer from normal cells, and is a critical step in the metastatic spread of the tumour. The invasion of the surrounding tissue involves the coordinated action between cancer cells, the ECM, the matrix degrading enzymes, and the epithelial-to-mesenchymal transition (EMT). This is a regulatory process through which epithelial cells (ECs) acquire mesenchymal characteristics and transform to mesenchymal-like cells (MCs). In this paper, we present a new mathematical model which describes the transition from a collective invasion strategy for the ECs to an individual invasion strategy for the MCs. We achieve this by formulating a coupled hybrid system consisting of partial and stochastic differential equations that describe the evolution of the ECs and the MCs, respectively. This approach allows one to reproduce in a very natural way fundamental qualitative features of the current biomedical understanding of cancer invasion that are not easily captured by classical modelling approaches, for example, the invasion of the ECM by self-generated gradients and the appearance of EC invasion islands outside of the main body of the tumour.

In this work an activator-depleted reaction-diffusion system is investigated on polar coordinates with the aim of exploring the relationship and the corresponding influence of domain size on the types of possible diffusion-driven instabilities. Quantitative relationships are found in the form of necessary conditions on the area of a disk-shape domain with respect to the diffusion and reaction rates for certain types of diffusion-driven instabilities to occur. Robust analytical methods are applied to find explicit expressions for the eigenvalues and eigenfunctions of the diffusion operator on a disk-shape domain with homogenous Neumann boundary conditions in polar coordinates. Spectral methods are applied using chebyshev non-periodic grid for the radial variable and Fourier periodic grid on the angular variable to verify the nodal lines and eigen-surfaces subject to the proposed analytical findings. The full classification of the parameter space in light of the bifurcation analysis is obtained and numerically verified by finding the solutions of the partitioning curves inducing such a classification. Furthermore, analytical results are found relating the area of (disk-shape) domain with reaction-diffusion rates in the form of necessary conditions for the different types of bifurcations. These results are on one hand presented in the form of mathematical theorems with rigorous proofs, and, on the other hand using finite element method, each claim of the corresponding theorems are verified by obtaining the theoretically predicted behaviour of the dynamics in the numerical simulations.

This paper explores the classification of parameter spaces for reaction-diffusion systems of two chemical species on stationary domains. The dynamics of the system are explored both in the absence and presence of diffusion. The parameter space is fully classified in terms of the types and stability of the uniform steady state. In the absence of diffusion the results on the classification of parameter space are supported by simulations of the corresponding vector-field and some trajectories around the uniform steady state. In the presence of diffusion, the main findings are the quantitative analysis relating the domain-size with the reaction and diffusion rates and their corresponding influence on the dynamics of the reaction-diffusion system when perturbed in the neighbourhood of the uniform steady state. Theoretical predictions are supported by numerical simulations both in the presence as well as in the absence of diffusion. Conditions on the domain size with respect to the diffusion and reaction rates are related to the types of diffusion-driven instabilities namely Turing, Hopf and Transcritical types of bifurcations. The first condition is an upper bound on the area of a rectangular domain in terms of the diffusion and reaction rates, which forbids the existence of Hopf and Transcritical types of bifurcations, yet allowing Turing type instability to occur. The second condition (necessary) is a lower bound on the domain size in terms of the reaction and diffusion rates to give rise to Hopf and Transcritical types of bifurcations as well as Turing instability.

We propose and analyse a novel surface finite element method that preserves the invariant regions of systems of semilinear parabolic equations on closed compact surfaces in $\mathbb{R}^3$ under discretisation. We also provide a fully-discrete scheme by applying the implicit-explicit (IMEX) Euler method in time. We prove the preservation of the invariant rectangles of the continuous problem under spatial and full discretizations. For scalar equations, these results reduce to the well-known discrete maximum principle. Furthermore, we prove optimal error bounds for the semi- and fully-discrete methods, that is the convergence rates are quadratic in the meshsize and linear in the timestep. Numerical experiments are provided to support the theoretical findings. In particular we provide examples in which, in the absence of lumping, the numerical solution violates the invariant region leading to blow-up due to the nature of the kinetics., Comment: 34 pages, 5 figures, 4 tables

We present a Bayesian methodology for infinite as well as finite dimensional parameter identification for partial differential equation models. The Bayesian framework provides a rigorous mathematical framework for incorporating prior knowledge on uncertainty in the observations and the parameters themselves, resulting in an approximation of the full probability distribution for the parameters, given the data. Although the numerical approximation of the full probability distribution is computationally expensive, parallelised algorithms can make many practically relevant problems computationally feasible. The probability distribution not only provides estimates for the values of the parameters, but also provides information about the inferability of parameters and the sensitivity of the model. This information is crucial when a mathematical model is used to study the outcome of real-world experiments. Keeping in mind the applicability of our approach to tackle real-world practical problems with data from experiments, in this initial proof of concept work, we apply this theoretical and computational framework to parameter identification for a well studied semilinear reaction-diffusion system with activator-depleted reaction kinetics, posed on evolving and stationary domains.

In this article we present a computational framework for isolating spatial patterns arising in the steady states of reaction-diffusion systems. Such systems have been used to model many different phenomena in areas such as developmental and cancer biology, cell motility and material science. Often one is interested in identifying parameters which will lead to a particular pattern. To attempt to answer this, we compute eigenpairs of the Laplacian on a variety of domains and use linear stability analysis to determine parameter values for the system that will lead to spatially inhomogeneous steady states whose patterns correspond to particular eigenfunctions. This method has previously been used on domains and surfaces where the eigenvalues and eigenfunctions are found analytically in closed form. Our contribution to this methodology is that we numerically compute eigenpairs on arbitrary domains and surfaces. Here we present various examples and demonstrate that mode isolation is straightforward especially for low eigenvalues. Additionally we see that if two or more eigenvalues are in a permissible range then the inhomogeneous steady state can be a linear combination of the respective eigenfunctions. Finally we show an example which suggests that pattern formation is robust on similar surfaces in cases that the surface either has or does not have a boundary.

Cell tracking algorithms which automate and systematise the analysis of time lapse image data sets of cells are an indispensable tool in the modelling and understanding of cellular phenomena. In this study we present a theoretical framework and an algorithm for whole cell tracking. Within this work we consider that "tracking" is equivalent to a dynamic reconstruction of the whole cell data (morphologies) from static image datasets. The novelty of our work is that the tracking algorithm is driven by a model for the motion of the cell. This model may be regarded as a simplification of a recently developed physically meaningful model for cell motility. The resulting problem is the optimal control of a geometric evolution law and we discuss the formulation and numerical approximation of the optimal control problem. The overall goal of this work is to design a framework for cell tracking within which the recovered data reflects the physics of the forward model. A number of numerical simulations are presented that illustrate the applicability of our approach.

Keratin are among the most abundant proteins in epithelial cells. Functions of the keratin network in cells are shaped by their dynamical organization. Using a collection of experimentally-driven mathematical models, different hypotheses for the turnover and transport of the keratin material in epithelial cells are tested. The interplay between turnover and transport and their effects on the keratin organization in cells are hence investigated by combining mathematical modeling and experimental data. Amongst the collection of mathematical models considered, a best model strongly supported by experimental data is identified. Fundamental to this approach is the fact that optimal parameter values associated with the best fit for each model are established. The best candidate among the best fits is characterized by the disassembly of the assembled keratin material in the perinuclear region and an active transport of the assembled keratin. Our study shows that an active transport of the assembled keratin is required to explain the experimentally observed keratin organization., Comment: 27 pages, 11 Figures

In this article we present robust, efficient and accurate fully implicit time-stepping schemes and nonlinear solvers for systems of reaction-diffusion equations. The applications of reaction-diffusion systems is abundant in the literature, from modelling pattern formation in developmental biology to cancer research, wound healing, tissue and bone regeneration and cell motility. Therefore, it is crucial that modellers, analysts and biologists are able to solve accurately and efficiently systems of highly nonlinear parabolic partial differential equations on complex stationary and sometimes continuously evolving domains and surfaces. The main contribution of our paper is the study of fully implicit schemes by use of the Newton method and the Picard iteration applied to the backward Euler, the Crank-Nicolson (and its modifications) and the fractional-step theta methods. Our results conclude that the fractional-step theta method coupled with a single Newton iteration at each timestep is as accurate as the fully adaptive Newton method; and both outperform the Picard iteration. In particular, the results strongly support the observation that a single Newton iteration is sufficient to yield as accurate results as those obtained by use of an adaptive Newton method. This is particularly advantageous when solving highly complex nonlinear partial differential equations on evolving domains and surfaces. To validate our theoretical results, various appropriate numerical experiments are exhibited on stationary planary domains and in the bulk of stationary surfaces., Comment: 19 pages, 8 figures

In this article we formulate new models for coupled systems of bulk-surface reaction-diffusion equations on stationary volumes. The bulk reaction-diffusion equations are coupled to the surface reaction-diffusion equations through linear Robin-type boundary conditions. We then state and prove the necessary conditions for diffusion-driven instability for the coupled system. Due to the nature of the coupling between bulk and surface dynamics, we are able to decouple the stability analysis of the bulk and surface dynamics. Under a suitable choice of model parameter values, the bulk reaction-diffusion system can induce patterning on the surface independent of whether the surface reaction-diffusion system produces or not, patterning. On the other hand, the surface reaction-diffusion system can not generate patterns everywhere in the bulk in the absence of patterning from the bulk reaction-diffusion system. For this case, patterns can only be induced in regions close to the surface membrane. Various numerical experiments are presented to support our theoretical findings. Our most revealing numerical result is that, Robin-type boundary conditions seem to introduce a boundary layer coupling the bulk and surface dynamics., Comment: 21 pages, 5 figures

We propose an adaptive finite element method to approximate the solutions to reaction-diffusion systems on time-dependent domains and surfaces. We derive a computable error estimator that provides an upper bound for the error in the semidiscrete (space) scheme. We reconcile our theoretical results with benchmark computations., Comment: 9 Pages, 4 Figures

We present and analyse an implicit-explicit timestepping procedure with finite element spatial approximation for a semilinear reaction-diffusion systems on evolving domains arising from biological models, such as Schnakenberg's (1979). We employ a Lagrangian formulation of the model equations which permits the error analysis for parabolic equations on a fixed domain but introduces technical difficulties, foremost the space-time dependent conductivity and diffusion. We prove optimal-order error estimates in the $\Lp{\infty}(0,T;\Lp{2}(\W))$ and $\Lp{2}(0,T;\Hil{1}(\W))$ norms, and a pointwise stability result. We remark that these apply to Eulerian solutions. Details on the implementation of the Lagrangian and the Eulerian scheme are provided. We also report on a numerical experiment for an application to pattern formation on an evolving domain.

We present global existence results for solutions of reaction-diffusion systems on evolving domains. Global existence results for a class of reaction-diffusion systems on fixed domains are extended to the same systems posed on spatially linear isotropically evolving domains. The results hold without any assumptions on the sign of the growth rate. The analysis is valid for many systems that commonly arise in the theory of pattern formation. We present numerical results illustrating our theoretical findings., Comment: 24 pages, 3 figures