Novel Aspects in Pattern Formation Arise from Coupling Turing Reaction-Diffusion and Chemotaxis.

Recent experimental studies on primary hair follicle formation and feather bud morphogenesis indicate a coupling between Turing-type diffusion driven instability and chemotactic patterning. Inspired by these findings we develop and analyse a mathematical model that couples chemotaxis to a reaction-diffusion system exhibiting diffusion-driven (Turing) instability. While both systems, reaction-diffusion systems and chemotaxis, can independently generate spatial patterns, we were interested in how the coupling impacts the stability of the system, parameter region for patterning, pattern geometry, as well as the dynamics of pattern formation. We conduct a classical linear stability analysis for different model structures, and confirm our results by numerical analysis of the system. Our results show that the coupling generally increases the robustness of the patterning process by enlarging the pattern region in the parameter space. Concerning time scale and pattern regularity, we find that an increase in the chemosensitivity can speed up the patterning process for parameters inside and outside of the Turing space, but generally reduces spatial regularity of the pattern. Interestingly, our analysis indicates that pattern formation can also occur when neither the Turing nor the chemotaxis system can independently generate pattern. On the other hand, for some parameter settings, the coupling of the two processes can extinguish the pattern formation, rather than reinforce it. These theoretical findings can be used to corroborate the biological findings on morphogenesis and guide future experimental studies. From a mathematical point of view, this work sheds a light on coupling classical pattern formation systems from the parameter space perspective.

A novel nonlocal partial differential equation model of endothelial progenitor cell cluster formation during the early stages of vasculogenesis.

The formation of new vascular networks is essential for tissue development and regeneration, in addition to playing a key role in pathological settings such as ischemia and tumour development. Experimental findings in the past two decades have led to the identification of a new mechanism of neovascularisation, known as cluster-based vasculogenesis, during which endothelial progenitor cells (EPCs) mobilised from the bone marrow are capable of bridging distant vascular beds in a variety of hypoxic settings in vivo. This process is characterised by the formation of EPC clusters during its early stages and, while much progress has been made in identifying various mechanisms underlying cluster formation, we are still far from a comprehensive description of such spatio-temporal dynamics. In order to achieve this, we propose a novel mathematical model of the early stages of cluster-based vasculogenesis, comprising of a system of nonlocal partial differential equations including key mechanisms such as endogenous chemotaxis, matrix degradation, cell proliferation and cell-to-cell adhesion. We conduct a linear stability analysis on the system and solve the equations numerically. We then conduct a parametric analysis of the numerical solutions of the one-dimensional problem to investigate the role of underlying dynamics on the speed of cluster formation and the size of clusters, measured via appropriate metrics for the cluster width and compactness. We verify the key results of the parametric analysis with simulations of the two-dimensional problem. Our results, which qualitatively compare with data from in vitro experiments, elucidate the complementary role played by endogenous chemotaxis and matrix degradation in the formation of clusters, suggesting chemotaxis is responsible for the cluster topology while matrix degradation is responsible for the speed of cluster formation. Our results also indicate that the nonlocal cell-to-cell adhesion term in our model, even though it initially causes cells to aggregate, is not sufficient to ensure clusters are stable over long time periods. Consequently, new modelling strategies for cell-to-cell adhesion are required to stabilise in silico clusters. We end the paper with a thorough discussion of promising, fruitful future modelling and experimental research perspectives.

Mechanical Models of Pattern and Form in Biological Tissues: The Role of Stress-Strain Constitutive Equations.

Mechanical and mechanochemical models of pattern formation in biological tissues have been used to study a variety of biomedical systems, particularly in developmental biology, and describe the physical interactions between cells and their local surroundings. These models in their original form consist of a balance equation for the cell density, a balance equation for the density of the extracellular matrix (ECM), and a force-balance equation describing the mechanical equilibrium of the cell-ECM system. Under the assumption that the cell-ECM system can be regarded as an isotropic linear viscoelastic material, the force-balance equation is often defined using the Kelvin-Voigt model of linear viscoelasticity to represent the stress-strain relation of the ECM. However, due to the multifaceted bio-physical nature of the ECM constituents, there are rheological aspects that cannot be effectively captured by this model and, therefore, depending on the pattern formation process and the type of biological tissue considered, other constitutive models of linear viscoelasticity may be better suited. In this paper, we systematically assess the pattern formation potential of different stress-strain constitutive equations for the ECM within a mechanical model of pattern formation in biological tissues. The results obtained through linear stability analysis and the dispersion relations derived therefrom support the idea that fluid-like constitutive models, such as the Maxwell model and the Jeffrey model, have a pattern formation potential much higher than solid-like models, such as the Kelvin-Voigt model and the standard linear solid model. This is confirmed by the results of numerical simulations, which demonstrate that, all else being equal, spatial patterns emerge in the case where the Maxwell model is used to represent the stress-strain relation of the ECM, while no patterns are observed when the Kelvin-Voigt model is employed. Our findings suggest that further empirical work is required to acquire detailed quantitative information on the mechanical properties of components of the ECM in different biological tissues in order to furnish mechanical and mechanochemical models of pattern formation with stress-strain constitutive equations for the ECM that provide a more faithful representation of the underlying tissue rheology.

A space-jump derivation for non-local models of cell-cell adhesion and non-local chemotaxis.

Cellular adhesion provides one of the fundamental forms of biological interaction between cells and their surroundings, yet the continuum modelling of cellular adhesion has remained mathematically challenging. In 2006, Armstrong et al. proposed a mathematical model in the form of an integro-partial differential equation. Although successful in applications, a derivation from an underlying stochastic random walk has remained elusive. In this work we develop a framework by which non-local models can be derived from a space-jump process. We show how the notions of motility and a cell polarization vector can be naturally included. With this derivation we are able to include microscopic biological properties into the model. We show that particular choices yield the original Armstrong model, while others lead to more general models, including a doubly non-local adhesion model and non-local chemotaxis models. Finally, we use random walk simulations to confirm that the corresponding continuum model represents the mean field behaviour of the stochastic random walk.

Structured models of cell migration incorporating molecular binding processes.

The dynamic interplay between collective cell movement and the various molecules involved in the accompanying cell signalling mechanisms plays a crucial role in many biological processes including normal tissue development and pathological scenarios such as wound healing and cancer. Information about the various structures embedded within these processes allows a detailed exploration of the binding of molecular species to cell-surface receptors within the evolving cell population. In this paper we establish a general spatio-temporal-structural framework that enables the description of molecular binding to cell membranes coupled with the cell population dynamics. We first provide a general theoretical description for this approach and then illustrate it with three examples arising from cancer invasion.

Robustness of a new molecular dynamics-finite element coupling approach for soft matter systems analyzed by uncertainty quantification.

Key parameters of a recently developed coarse-grained molecular dynamics-finite element coupling approach have been analyzed in the framework of uncertainty quantification (UQ). We have employed a polystyrene sample for the case study. The new hybrid approach contains several parameters which cannot be determined on the basis of simple physical arguments. Among others, this includes the so-called anchor points as information transmitters between the particle-based molecular dynamics (MD) domain and the surrounding finite element continuum, the force constant between polymer beads and anchor points, the number of anchor points, and the relative sizes of the MD core domain and the surrounding dissipative particle dynamics domain. Polymer properties such as density, radius of gyration, end-to-end distance, and radial distribution functions are calculated as a function of the above model parameters. The influence of these input parameters on the resulting polymer properties is studied by UQ. Our analysis shows that the hybrid method is highly robust. The variation of polymer properties of interest as a function of the input parameters is weak.

Mathematical modelling of cancer invasion: implications of cell adhesion variability for tumour infiltrative growth patterns.

Cancer invasion, recognised as one of the hallmarks of cancer, is a complex, multiscale phenomenon involving many inter-related genetic, biochemical, cellular and tissue processes at different spatial and temporal scales. Central to invasion is the ability of cancer cells to alter and degrade an extracellular matrix. Combined with abnormal excessive proliferation and migration which is enabled and enhanced by altered cell-cell and cell-matrix adhesion, the cancerous mass can invade the neighbouring tissue. Along with tumour-induced angiogenesis, invasion is a key component of metastatic spread, ultimately leading to the formation of secondary tumours in other parts of the host body. In this paper we explore the spatio-temporal dynamics of a model of cancer invasion, where cell-cell and cell-matrix adhesion is accounted for through non-local interaction terms in a system of partial integro-differential equations. The change of adhesion properties during cancer growth and development is investigated here through time-dependent adhesion characteristics within the cell population as well as those between the cells and the components of the extracellular matrix. Our computational simulation results demonstrate a range of heterogeneous dynamics which are qualitatively similar to the invasive growth patterns observed in a number of different types of cancer, such as tumour infiltrative growth patterns (INF).

Mathematical modeling of cancer cell invasion of tissue: biological insight from mathematical analysis and computational simulation.

The ability of cancer cells to break out of tissue compartments and invade locally gives solid tumours a defining deadly characteristic. One of the first steps of invasion is the remodelling of the surrounding tissue or extracellular matrix (ECM) and a major part of this process is the over-expression of proteolytic enzymes, such as the urokinase-type plasminogen activator (uPA) and matrix metalloproteinases (MMPs), by the cancer cells to break down ECM proteins. Degradation of the matrix enables the cancer cells to migrate through the tissue and subsequently to spread to secondary sites in the body, a process known as metastasis. In this paper we undertake an analysis of a mathematical model of cancer cell invasion of tissue, or ECM, which focuses on the role of the urokinase plasminogen activation system. The model consists of a system of five reaction-diffusion-taxis partial differential equations describing the interactions between cancer cells, uPA, uPA inhibitors, plasmin and the host tissue. Cancer cells react chemotactically and haptotactically to the spatio-temporal effects of the uPA system. The results obtained from computational simulations carried out on the model equations produce dynamic heterogeneous spatio-temporal solutions and using linear stability analysis we show that this is caused by a taxis-driven instability of a spatially homogeneous steady-state. Finally we consider the biological implications of the model results, draw parallels with clinical samples and laboratory based models of cancer cell invasion using three-dimensional invasion assay, and go on to discuss future development of the model.

Uncertainty Quantification Guided Parameter Selection in a Fully Coupled Molecular Dynamics-Finite Element Model of the Mechanical Behavior of Polymers.

The objective of investigating macroscopic polymer properties with a low computing cost and a high resolution has led to the development of efficient hybrid simulation tools. Systems generated from such simulation tools can fail in service if the effect of uncertainty of model inputs on its outputs is not accounted for. This work focuses on quantifying the effect of parametric uncertainty in our coarse-grained molecular dynamics-finite element coupling approach using uncertainty quantification. We consider uniaxial deformation simulations of a polystyrene sample at T = 100 K in our study. Parametric uncertainty is assumed to originate from parameters in the molecular dynamics model with a nonperiodic boundary (the force constant between polymer beads and anchor points, the number of anchor points, and the size of the surrounding dissipative particle dynamics domain) and a parameter to blend the energies of particles and continuum (weighting factor). Key issues that arise in uncertainty quantification are discussed on the basis of the quantities of interest including mass density, end-to-end distance, and radial distribution function. This work reveals the influence of key input parameters on the properties of polymer structure and facilitates the determination of those parameters in the application of this hybrid molecular dynamics-finite element approach.

Mathematical modeling in wound healing, bone regeneration and tissue engineering.

The processes of wound healing and bone regeneration and problems in tissue engineering have been an active area for mathematical modeling in the last decade. Here we review a selection of recent models which aim at deriving strategies for improved healing. In wound healing, the models have particularly focused on the inflammatory response in order to improve the healing of chronic wound. For bone regeneration, the mathematical models have been applied to design optimal and new treatment strategies for normal and specific cases of impaired fracture healing. For the field of tissue engineering, we focus on mathematical models that analyze the interplay between cells and their biochemical cues within the scaffold to ensure optimal nutrient transport and maximal tissue production. Finally, we briefly comment on numerical issues arising from simulations of these mathematical models.

Angiogenesis in bone fracture healing: a bioregulatory model.

The process of fracture healing involves the action and interaction of many cells, regulated by biochemical and mechanical signals. Vital to a successful healing process is the restoration of a good vascular network. In this paper, a continuous mathematical model is presented that describes the different fracture healing stages and their response to biochemical stimuli only (a bioregulatory model); mechanoregulatory effects are excluded here. The model consists of a system of nonlinear partial differential equations describing the spatiotemporal evolution of concentrations and densities of the cell types, extracellular matrix types and growth factors indispensable to the healing process. The model starts after the inflammation phase, when the fracture callus has already been formed. Cell migration is described using not only haptokinetic, but also chemotactic and haptotactic influences. Cell differentiation is controlled by the presence of growth factors and sufficient vascularisation. Matrix synthesis and growth factor production are controlled by the local cell and matrix densities and by the local growth factor concentrations. Numerical simulations of the system, using parameter values based on experimental data obtained from literature, are presented. The simulation results are corroborated by comparison with experimental data from a standardised rodent fracture model. The results of sensitivity analyses on the parameter values as well as on the boundary and initial conditions are discussed. Numerical simulations of compromised healing situations showed that the establishment of a vascular network in response to angiogenic growth factors is a key factor in the healing process. Furthermore, a correct description of cell migration is also shown to be essential to the prediction of realistic spatiotemporal tissue distribution patterns in the fracture callus. The mathematical framework presented in this paper can be an important tool in furthering the understanding of the mechanisms causing compromised healing and can be applied in the design of future fracture healing experiments.

Derivation of the mesoscopic elasticity tensor of cortical bone from quantitative impedance images at the micron scale.

This paper addresses the relationships between the microscopic properties of bone and its elasticity at the millimetre scale, or mesoscale. A method is proposed to estimate the mesoscale properties of cortical bone based on a spatial distribution of acoustic properties at the microscopic scale obtained with scanning acoustic microscopy. The procedure to compute the mesoscopic stiffness tensor involves (i) the segmentation of the pores to obtain a realistic model of the porosity; (ii) the construction of a field of anisotropic elastic coefficients at the microscopic scale which reflects the heterogeneity of the bone matrix; (iii) finite element computations of mesoscopic homogenized properties. The computed mesoscopic properties compare well with available experimental data. It appears that the tissue anisotropy at the microscopic level has a major effect on the mesoscopic anisotropy and that assuming the pores filled with an incompressible fluid or, alternatively, empty, leads to significantly different mesoscopic properties.

Mathematical modeling of fracture healing in mice: comparison between experimental data and numerical simulation results.

The combined use of experimental and mathematical models can lead to a better understanding of fracture healing. In this study, a mathematical model, which was originally established by Bailón-Plaza and van der Meulen (J Theor Biol 212:191-209, 2001), was applied to an experimental model of a semi-stabilized murine tibial fracture. The mathematical model was implemented in a custom finite volumes code, specialized in dealing with the model's requirements of mass conservation and non-negativity of the variables. A qualitative agreement between the experimentally measured and numerically simulated evolution in the cartilage and bone content was observed. Additionally, an extensive parametric study was conducted to assess the influence of the model parameters on the simulation outcome. Finally, a case of pathological fracture healing and its treatment by administration of growth factors was modeled to demonstrate the potential therapeutic value of this mathematical model.

To what extent can cortical bone millimeter-scale elasticity be predicted by a two-phase composite model with variable porosity?

An evidence gap exists in fully understanding and reliably modeling the variations in elastic anisotropy that are observed at the millimeter scale in human cortical bone. The porosity (pore volume fraction) is known to account for a large part, but not all, of the elasticity variations. This effect may be modeled by a two-phase micromechanical model consisting of a homogeneous matrix pervaded by cylindrical pores. Although this model has been widely used, it lacks experimental validation. The aim of the present work is to revisit experimental data (elastic coefficients, porosity) previously obtained from 21 cortical bone specimens from the femoral mid-diaphysis of 10 donors and test the validity of the model by proposing a detailed discussion of its hypotheses. This includes investigating to what extent the experimental uncertainties, pore network modeling, and matrix elastic properties influence the model's predictions. The results support the validity of the two-phase model of cortical bone which assumes that the essential source of variations of elastic properties at the millimeter-scale is the volume fraction of vascular porosity. We propose that the bulk of the remaining discrepancies between predicted stiffness coefficients and experimental data (RMSE between 6% and 9%) is in part due to experimental errors and part due to small variations of the extravascular matrix properties. More significantly, although most of the models that have been proposed for cortical bone were based on several homogenization steps and a large number of variable parameters, we show that a model with a single parameter, namely the volume fraction of vascular porosity, is a suitable representation for cortical bone. The results could provide a guide to build specimen-specific cortical bone models. This will be of interest to analyze the structure-function relationship in bone and to design bone-mimicking materials.

On the elastic properties of mineralized turkey leg tendon tissue: multiscale model and experiment.

The key parameters influencing the elastic properties of the mineralized turkey leg tendon (MTLT) were investigated. Two structurally different tissue types appearing in the MTLT were considered: circumferential and interstitial tissue. These differ in their amount of micropores and their average diameter of the mineralized collagen fibril bundles. A multiscale model representing the apparent elastic stiffness tensor of MTLT tissue was developed using the Mori-Tanaka and the self-consistent homogenization schemes. The volume fraction of mineral (hydroxyapatite) in the fibril bundle, [Formula: see text], and the tissue microporosity are the variables of the model. The MTLT model was analyzed performing a global sensitivity analysis (Elementary Effects method) and a parametric study. The stiffnesses parallel (axial) and perpendicular (transverse) to the MTLT long axis were the only significantly sensitive components of the apparent stiffness tensor of MTLT tissue. The most important parameters influencing these apparent stiffnesses are [Formula: see text], tissue microporosity, as well as shape and distribution of the minerals in the fibril bundle (intra- vs. interfibrillar). The predicted apparent stiffness was converted to acoustic impedance for model validation. From measurements on embedded MTLT samples, including 50- and 200-MHz scanning acoustic microscopy as well as synchrotron radiation micro-computed tomography, we obtained site-matched acoustic impedance and [Formula: see text] data of circumferential and interstitial tissue. The experimental and the model data compare very well for both tissue types (relative error 6-8 %).

Spatial distribution of tissue level properties in a human femoral cortical bone.

The mechanical properties of cortical bone are determined by a combination bone tissue composition, and structure at several hierarchical length scales. In this study the spatial distribution of tissue level properties within a human femoral shaft has been investigated. Cylindrically shaped samples (diameter: 4.4mm, N=56) were prepared from cortical regions along the entire length (20-85% of the total femur length), and around the periphery (anterior, medial, posterior and lateral quadrants). The samples were analyzed using scanning acoustic microscopy (SAM) at 50MHz and synchrotron radiation micro computed tomography (SRµCT). For all samples the average cortical porosity (Ct.Po), tissue elastic coefficients (c(ij)) and the average tissue degree of mineralization (DMB) were determined. The smallest coefficient of variation was observed for DMB (1.8%), followed by BV/TV (5.4%), c(ij) (8.2-45.5%), and Ct.Po (47.5%). Different variations with respect to the anatomical position were found for DMB, Ct.Po and c(ij). These data address the anatomical variations in anisotropic elastic properties and link them to tissue mineralization and porosity, which are important input parameters for numerical multi-scale bone models.

A determination of the minimum sizes of representative volume elements for the prediction of cortical bone elastic properties.

At its highest level of microstructural organization-the mesoscale or millimeter scale-cortical bone exhibits a heterogeneous distribution of pores (Haversian canals, resorption cavities). Multi-scale mechanical models rely on the definition of a representative volume element (RVE). Analytical homogenization techniques are usually based on an idealized RVE microstructure, while finite element homogenization using high-resolution images is based on a realistic RVE of finite size. The objective of this paper was to quantify the size and content of possible cortical bone mesoscale RVEs. RVE size was defined as the minimum size: (1) for which the apparent (homogenized) stiffness tensor becomes independent of the applied boundary conditions or (2) for which the variance of elastic properties for a set of microstructure realizations is sufficiently small. The field of elastic coefficients and microstructure in RVEs was derived from one acoustic microscopy image of a human femur cortical bone sample with an overall porosity of 8.5%. The homogenized properties of RVEs were computed with a finite element technique. It was found that the size of the RVE representative of the overall tissue is about 1.5 mm. Smaller RVEs (~0.5 mm) can also be considered to estimate local mesoscopic properties that strongly depend on the local pores volume fraction. This result provides a sound basis for the application of homogenization techniques to model the heterogeneity of cortical microstructures. An application of the findings to estimate elastic properties in the case of a porosity gradient is briefly presented.

A hybrid bioregulatory model of angiogenesis during bone fracture healing.

Bone fracture healing is a complex process in which angiogenesis or the development of a blood vessel network plays a crucial role. In this paper, a mathematical model is presented that simulates the biological aspects of fracture healing including the formation of individual blood vessels. The model consists of partial differential equations, several of which describe the evolution in density of the most important cell types, growth factors, tissues and nutrients. The other equations determine the growth of blood vessels as a result of the movement of leading endothelial (tip) cells. Branching and anastomoses are accounted for in the model. The model is applied to a normal fracture healing case and subjected to a sensitivity analysis. The spatiotemporal evolution of soft tissues and bone, as well as the development of a blood vessel network are corroborated by comparison with experimental data. Moreover, this study shows that the proposed mathematical framework can be a useful tool in the research of impaired healing and the design of treatment strategies.