Axolotl salamanders regrow entire limbs throughout life using a molecular machinery similar to that used in the embryonic development of our limbs. A series of genetically-regulated molecular markers, known as morphogens, coordinate key mechanisms of this process. But exactly how morphogen patterning drives skeletal joint formation is not fully understood yet. Computational models provide a means to explore the factors influencing this process. In addition, we have access to three-dimensional (3D) microscopy images of morphogen expression in regenerating axolotl forelimbs obtained using a novel technique [1], which can help inform such models. But to effectively abstract key concepts and test hypotheses about how biological processes work using computational models, we must have a thorough grasp on how model parameters and conditions influence predicted outcomes.

One of the prevailing theories to explain the formation of morphogen patterns was proposed by the mathematician Alan Turing [2]. Starting from a nearly uniform initial state, and through suitable nonlinear interactions between reacting and diffusing morphogens, stable spatial patterns, known as Turing patterns, are obtained. Turing patterns have been used to investigate the formation of skeletal limb structures, but previous studies only used one or two-dimensional models, which have limited applicability due to the joint’s asymmetrical structure. Extending to 3D is challenging due to the wider variety of patterns and increased complexity.

We use linear stability analysis and finite element modeling to predict Turing pattern emergence in a 3D generic domain using the Schnakenberg activator-substrate model [3]. We provide a framework to identify the critical factors necessary for specific morphogen pattern emergence, and explore the role of initial conditions, model parameters and domain size on the predictions. We have observed that initial conditions on the activator have a stronger impact on the final pattern than initial conditions on the substrate, which has important implications for modeling morphogen expression in joint formation.

Our findings establish the groundwork for upcoming research, where we aim to employ the model developed here to predict the morphogen expression resulting from the axolotl joint formation experiments. Only through a comprehensive understanding of the factors influencing pattern emergence will we be able to successfully use Turing patterns as a model for morphogen patterning in joint formation.

References

[1] Lovely AM et al. HCR-FISH in Ambystoma mexicanum Tissue. In: Salamanders Methods in Molecular Biology. New York, 2023. p.109–22, doi:10.1007/978-1-0716-2659-7_6.

[2] Turing A. The chemical basis of morphogenesis. Philos Trans R Soc London, 1952; 237(641):37-72, doi:10.1007/BF02459572.

[3] Ben Tahar S at al. Turing pattern prediction in three-dimensional domains: the role of initial conditions and growth. Preprint on bioRxiv, 2023; doi: 10.1101/2023.03.29.534782.

Cells are constantly interacting with their environment, adapting their behavior in response to different stimuli and environmental conditions. This cellular adaptation occurs via changes in gene expression derived from changes in the physiological environment, giving rise to phenotypic plasticity. Phenotypic plasticity plays a key role in many steps of cancer progression, to the point that it has been included among the Hallmarks of Cancer. Indeed, it partly explains some of the most characteristic features of cancer, such as metastasis or drug resistance, probably the main challenges for improving cancer prognosis. Hence, understanding the mechanisms that trigger cellular adaptation is crucial for advancing in the study of this disease and its eventual treatments

In this line, mathematical models and simulation have great potential for gaining insight into complex cell process such as adaptation, and testing hypotheses regarding the effects of different environmental conditions in the adaptive response of tumors. In this work, we present a novel mathematical framework to model cellular adaptation and phenotypic heterogeneity in cell populations interacting with their environment. The model is based on the concept of internal variables, which are used to model cell state and regulate cell behaviour, in addition to Partial Differential Equations (PDEs) to describe the evolution of cell populations and the spatial concentration of chemical species from a continuum point of view. The proposed approach allows to consider not only cell response to environmental changes, but also reversibility and inheritance typical of phenotypic changes.

After presentation of the model, we particularize it to the case of glioblastoma (GBM) adaptation to hypoxia. GBM is the most common and lethal primary brain tumor. It has a dismal prognosis with a 5-year survival rate of 5%. The poor response of GBM to treatment is a consequence of its intrinsic and acquired drug resistance. This resistance is enhanced by hypoxia, a defining feature of GBM. Therefore, it is important to study how hypoxia governs cellular adaptation in GBM to improve our understanding on treatment response and, eventually, GBM prognosis.

The objective of this study is analyzing the potential of the derived model for capturing important biological trends in GBM evolution and adaptation under variable oxygen concentrations. After an extensive parametric analysis, different oxygenation conditions are tested. The results show the flexibility of the model for capturing the variability existing among GBM tumors. The model is also able to capture some observed experimental trends, such as the increased aggressiveness and resilience of tumors undergoing cyclic hypoxia.

In short, the developed framework presents an alternative to modelling cellular adaptation and may, with suitable validation, help in designing predictive tools as well as in silico clinical trials.

Last few years have brought into broader scientific audience examples of very robust and prominent periodic spatio-temporal patterns in biological systems [1,2]. Until recently, patterns such as spirals and planar wave trains were observed almost exclusively in pure chemical nonlinear systems. Reaction-diffusion equations are commonly used for understanding such phenomena. In these systems, onset of periodic wave trains is described by the so called wave instability (WI). Here we present generic three component activator-depleted substrate-inhibitor model which is a minimal model for patterns created by small GTPase signalling molecule RhoA. The model reproduces patterns which are observed experimentally [1] and explains some quantitative relationships between amplitude, temporal period and concentrations of some components.

Activity of RhoA may also induce contraction of cell surface (cortex) by activation of a motor protein myosin, which may result in advection of chemical components. Hence, our system becomes a reaction-diffusion-advection system with temporally and spatially dependent velocity field. We are interested in modification of chemical instabilities by biomechanical modes and interactions of those modes in simple reaction-diffusion-active gel framework. In this approach cell cortex is treated as a visco-elastic fluid in which active stress is controlled locally by the chemical subsystem [3]. We found that biomechanical coupling may increase parameter range in which spatially homogenous solution is unstable and periodic spatial patterns are possible. We show that interactions between chemical and biomechanical modes may lead to the change of spatiotemporal characteristics of patterns created by the chemical subsystem, their destabilization as well as formation of new patterns. Moreover, we discuss the role of the delay between RhoA activity and active stress. 2D simulations on the plane and stability analysis of 1D homogenous and nonhomogeneous solutions are presented.

[1] A. Michaud, M. Leda, Z. Swider, S. Kim, J. He, J. Valley, J. Huisken, J. Landino, A. Goryachev, G. von Dassow, W. Bement, A versatile cortical pattern-forming circuit based on Rho, F-actin, Ect2, and RGA-3/4. J. Cell Biol., 221(8):e202203017, (2022)

[2] W. Bement, M. Leda, A. Moe, A. Kita, M. Larson, A. Golding, C. Pfeuti, K-C. Su, A. Miller, A. Goryachev, G. von Dassow, Activator-inhibitor coupling between Rho signalling and actin assembly makes the cell cortex an excitable medium, Nature Cell Biol., 17(11), 1471 – 1483 (2015).

[3] Frank Jülicher et al, Hydrodynamic theory of active matter, Rep. Prog. Phys. 81, 076601 (2018).

Red blood cells (RBCs) are known as important formed elements in the blood. The functions of RBCs are to transport oxygen from the lungs to the tissues and metabolic waste, carbon dioxide, from the tissues to the lungs and to maintain systemic acid/base equilibrium. Along with that, RBC is known to release adenosine triphosphate (ATP) when it is subjected to shear stress. Subsequently, ATP molecules bind to purinergic receptors, to activate a cascade of biochemical reactions in endothelial cells (ECs) to release sequestrated ubiquitous calcium ion from endoplasmic reticulum (ER). The response of EC to ATP is in the form of transient calcium. Calcium is well known to regulate the activity of many enzymes in order to maintain the cellular homeostasis. Nevertheless, unregulated free calcium concentration could lead to serious pathological conditions such as cell death. In order to avoid such consequences, EC manages to maintain its physiological calcium concentration in the presence of ATP using ATPdriven pumps present in the plasma/ER membrane and the desensitization of the purinergic receptors. In order to understand how ATP released from RBCs affect the intracellular homeostasis in a vascular wall, we firstly developed a minimal calcium model, which guarantees the intracellular calcium homeostasis. Secondly, we couple it to RBC flow in a two-dimensional channel for a given imposed parabolic flow. In simulation, we use immersed boundary-lattice Boltzmann method (IB-LBM) to solve the fluid flow and the ATP release from RBC. We carried out several simulations varying flow strength, channel width, and concentrations of RBC (hematocrit) in order to emulate the blood flow in microcirculation. The endothelium helps maintaining a steady ATP concentration, avoiding the abnormal rise in the ATP concentration released from RBCs. With varying flow strength and hematocrit for a given channel width, we found that the ATP concentration and the cytoplasmic transient concentration increase with increase in the flow strength as well as hematocrit, and this leads to the cytoplasmic transient time decrease. Due to the relatively small peak times and amplitudes of cytoplasmic calcium at high flow strength as compared to that at low flow strength for all hematocrit, there is a possibility of calcium propagation from the high flow strength region to the low flow strength region for the coordination of cellular functions. Similarly, we observed a possibility of calcium propagation from low confined channels to the medium or highly confined channels for all hematocrit for a given flow strength. It would be interesting to carry out a further study in a vascular network in order to get more insights on the calcium propagation as each branch of vascular network may have unequal concentrations of RBC and the flow strength.