A fast computational model for circulatory dynamics: effects of left ventricle–aorta coupling

The course of diseases such as hypertension, systolic heart failure and heart failure with a preserved ejection fraction is affected by interactions between the left ventricle (LV) and the vasculature. To study these interactions, a computationally efficient, biophysically based mathematical model for the circulatory system is presented. In a four-chamber model of the heart, the LV is represented by a previously described low-order, wall volume-preserving model that includes torsion and base-to-apex and circumferential wall shortening and lengthening, and the other chambers are represented using spherical geometries. Active and passive myocardial mechanics of all four chambers are included. The cardiac model is coupled with a wave propagation model for the aorta and a closed lumped-parameter circulation model. Parameters for the normal heart and aorta are determined by fitting to experimental data. Changes in the timing and magnitude of pulse wave reflections by the aorta are demonstrated with changes in compliance and taper of the aorta as seen in aging (decreased compliance, increased diameter and length), and resulting effects on LV pressure–volume loops and LV fiber stress and sarcomere shortening are predicted. Effects of aging of the aorta combined with reduced LV contractile force (failing heart) are examined. In the failing heart, changes in aortic properties with aging affect stroke volume and sarcomere shortening without appreciable augmentation of aortic pressure, and the reflected pressure wave contributes an increased proportion of aortic pressure.

     

Comparative chondrogenesis of human cells in a 3D integrated experimental–computational mechanobiology model

We present an integrated experimental–computational mechanobiology model of chondrogenesis. The response of human articular chondrocytes to culture medium perfusion, versus perfusion associated with cyclic pressurisation, versus non-perfused culture, was compared in a pellet culture model, and multiphysic computation was used to quantify oxygen transport and flow dynamics in the various culture conditions. At 2 weeks of culture, the measured cell metabolic activity and the matrix content in collagen type II and aggrecan were greatest in the perfused+pressurised pellets. The main effects of perfusion alone, relative to static controls, were to suppress collagen type I and GAG contents, which were greatest in the non-perfused pellets. All pellets showed a peripheral layer of proliferating cells, which was thickest in the perfused pellets, and most pellets showed internal gradients in cell density and matrix composition. In perfused pellets, the computed lowest oxygen concentration was 0.075 mM (7.5% tension), the maximal oxygen flux was 477.5 nmol/m²/s and the maximal fluid shear stress, acting on the pellet surface, was 1.8 mPa (0.018 dyn/cm²). In the non-perfused pellets, the lowest oxygen concentration was 0.003 mM (0.3% tension) and the maximal oxygen flux was 102.4 nmol/m²/s. A local correlation was observed, between the gradients in pellet properties obtained from histology, and the oxygen fields calculated with multiphysic simulation. Our results show up-regulation of hyaline matrix protein production by human chondrocytes in response to perfusion associated with cyclic pressurisation. These results could be favourably exploited in tissue engineering applications.     

A computational reaction–diffusion model for biosynthesis and linking of cartilage extracellular matrix in cell-seeded scaffolds with varying porosity

Biomechanics and Modeling in Mechanobiology - Cartilage tissue engineering is commonly initiated by seeding cells in porous materials such as hydrogels or scaffolds. Under optimal conditions, the...     

Emergence of orientation-selective inhibition in the primary visual cortex: a Bayes–Markov computational model

The recent consensus is that virtually all aspects of response selectivity exhibited by the primary visual cortex are either created or sharpened by cortical inhibitory interneurons. Experimental studies have shown that there are cortical inhibitory cells that are driven by geniculate cells and that, like their cortical excitatory counterparts, are orientation selective, though less sharply tuned. The main goal of this article is to demonstrate how orientation-selective inhibition might be created by the circuitry of the primary visual cortex (striate cortex, V1) from its nonoriented geniculate inputs. To fulfill this goal, first, a Bayes–Markov computational model is developed for the V1 area dedicated to foveal vision. The developed model consists of three parts: (i) a two-layered hierarchical Markov random field that is assumed to generate the activity patterns of the geniculate and cortical inhibitory cells, (ii) a Bayesian computational goal that is formulated based on the maximum a posteriori (MAP) estimation principle, and (iii) an iterative, deterministic, parallel algorithm that leads the cortical circuitry to achieve its assigned computational goal. The developed model is not fully LGN driven and it is not implementable by the neural machinery of V1. The model, then, is transformed into a fully LGN-driven and physiologically plausible form. Computer simulation is used to demonstrate the performance of the developed models.     

A computational framework for modeling cell–matrix interactions in soft biological tissues

Biomechanics and Modeling in Mechanobiology - Living soft tissues appear to promote the development and maintenance of a preferred mechanical state within a defined tolerance around a so-called set...     

A computational analysis of ordering in ABCHN at the nematic–isotropic transition

A computational analysis of ordering in the nematogenic compound 4-alkenyl bicyclohexylnitrile has been carried out based on quantum mechanics and intermolecular forces. The evaluation of atomic charge and dipole moment at each atomic center has been carried out using the complete neglect differential overlap (CNDO/2) method. Modified Rayleigh–Schrodinger perturbation theory along with a multicentered- multipole expansion method has been employed to evaluate long-range intermolecular interactions, while a ‘6-exp‘ potential function has been assumed for short-range interactions. The total interaction energy values obtained through these computations were used to calculate the probability of each configuration at room temperature (300 K), the nematic–isotropic transition temperature (364.7 K) and above transition temperature (450 K) using the Maxwell–Boltzmann formula. The various possible configurations during the different modes (i.e., stacking, in-plane and terminal) of interactions have been studied in terms of variation of probability due to small departures from the most probable configurations. An attempt has been made to analyze the characteristic features of liquid crystallinity in terms of their relative order with molecular parameters introduced in this paper.     

A new method for choosing the computational cell in stochastic reaction–diffusion systems

How to choose the computational compartment or cell size for the stochastic simulation of a reaction–diffusion system is still an open problem, and a number of criteria have been suggested. A generalized measure of the noise for finite-dimensional systems based on the largest eigenvalue of the covariance matrix of the number of molecules of all species has been suggested as a measure of the overall fluctuations in a multivariate system, and we apply it here to a discretized reaction–diffusion system. We show that for a broad class of first-order reaction networks this measure converges to the square root of the reciprocal of the smallest mean species number in a compartment at the steady state. We show that a suitably re-normalized measure stabilizes as the volume of a cell approaches zero, which leads to a criterion for the maximum volume of the compartments in a computational grid. We then derive a new criterion based on the sensitivity of the entire network, not just of the fastest step, that predicts a grid size that assures that the concentrations of all species converge to a spatially-uniform solution. This criterion applies for all orders of reactions and for reaction rate functions derived from singular perturbation or other reduction methods, and encompasses both diffusing and non-diffusing species. We show that this predicts the maximal allowable volume found in a linear problem, and we illustrate our results with an example motivated by anterior-posterior pattern formation in Drosophila, and with several other examples.     

A general reaction–diffusion model of acidity in cancer invasion

We model the metabolism and behaviour of a developing cancer tumour in the context of its microenvironment, with the aim of elucidating the consequences of altered energy metabolism. Of particular interest is the Warburg Effect, a widespread preference in tumours for cytosolic glycolysis rather than oxidative phosphorylation for glucose breakdown, as yet incompletely understood. We examine a candidate explanation for the prevalence of the Warburg Effect in tumours, the acid-mediated invasion hypothesis, by generalising a canonical non-linear reaction–diffusion model of acid-mediated tumour invasion to consider additional biological features of potential importance. We apply both numerical methods and a non-standard asymptotic analysis in a travelling wave framework to obtain an explicit understanding of the range of tumour behaviours produced by the model and how fundamental parameters govern the speed and shape of invading tumour waves. Comparison with conclusions drawn under the original system—a special case of our generalised system—allows us to comment on the structural stability and predictive power of the modelling framework.     

Replicability or reproducibility? On the replication crisis in computational neuroscience and sharing only relevant detail

Replicability and reproducibility of computational models has been somewhat understudied by “the replication movement.” In this paper, we draw on methodological studies into the replicability of psychological experiments and on the mechanistic account of explanation to analyze the functions of model replications and model reproductions in computational neuroscience. We contend that model replicability, or independent researchers' ability to obtain the same output using original code and data, and model reproducibility, or independent researchers' ability to recreate a model without original code, serve different functions and fail for different reasons. This means that measures designed to improve model replicability may not enhance (and, in some cases, may actually damage) model reproducibility. We claim that although both are undesirable, low model reproducibility poses more of a threat to long-term scientific progress than low model replicability. In our opinion, low model reproducibility stems mostly from authors' omitting to provide crucial information in scientific papers and we stress that sharing all computer code and data is not a solution. Reports of computational studies should remain selective and include all and only relevant bits of code.     

Elasticity and geometry: a computational model of the Heineke–Mikulicz strictureplasty

Crohn’s disease is a challenging inflammatory process with a propensity for focal gastro-intestinal tract inflammation and stricture. Surgically, Crohn’s is often treated with resection. However, a subtype of diffuse disease with multiple strictures is treated by strictureplasty procedures in hope of avoiding short-gut syndrome. Prior work by Pocivavsek et al. defined the geometry of a Heineke–Mikulicz strictureplasty. Here, we bring this analysis one step closer to clinical and biological relevance by calculating the mechanical stresses and strains that the strictureplasty deformation generates on a model intestinal wall. The small bowel is simulated as a linearly elastic isotropic deformable cylindrical shell using finite element modeling. Data show a divergence in elastic response between the anti-mesenteric and mesenteric halves. The anti-mesenteric surface shows a bending dominated elastic response that correlates with the prior purely geometric analysis. However, the mesenteric half is not a neutral bystander during strictureplasty formation, as geometric arguments predict. Strong in-plane stretching strains develop in a rim around the image of the transverse closure, which may impact local perfusion and serve as sites of disease recurrence. Lastly, nearly all the deformation energy is stored in the central vertex stitch, placing this part at highest risk of dehiscence. This study enhances our understanding of mechanical response in complex nonlinear cylindrical geometries like the surgically manipulated intestinal tract. The developed framework serves as a platform for future addition of more complex clinically relevant parameters to our model, including real tissue properties, anisotropy, blood supply modeling, and patient deriver anatomic factors.     

A model of a MAPK•substrate complex in an active conformation: a computational and experimental approach

The mechanisms by which MAP kinases recognize and phosphorylate substrates are not completely understood. Efforts to understand the mechanisms have been compromised by the lack of MAPK-substrate structures. While MAPK-substrate docking is well established as a viable mechanism for bringing MAPKs and substrates into close proximity the molecular details of how such docking promotes phosphorylation is an unresolved issue. In the present study computer modeling approaches, with restraints derived from experimentally known interactions, were used to predict how the N-terminus of Ets-1 associates with ERK2. Interestingly, the N-terminus does not contain a consensus-docking site ((R/K)_2-3-X_2-6-Φ_A-X-Φ_B, where Φ is aliphatic hydrophobic) for ERK2. The modeling predicts that the N-terminus of Ets-1 makes important contributions to the stabilization of the complex, but remains largely disordered. The computer-generated model was used to guide mutagenesis experiments, which support the notion that Leu-11 and possibly Ile-13 and Ile-14 of Ets-1 1-138 (Ets) make contributions through binding to the hydrophobic groove of the ERK2 D-recruiting site (DRS). Based on the modeling, a consensus-docking site was introduced through the introduction of an arginine at residue 7, to give the consensus ⁷RK-X₂-Φ_A-X-Φ_B ¹³. This results in a 2-fold increase in k _cat/K _m for the phosphorylation of Ets by ERK2. Similarly, the substitution of the N-terminus for two different consensus docking sites derived from Elk-1 and MKK1 also improves k _cat/K _m by two-fold compared to Ets. Disruption of the N-terminal docking through deletion of residues 1-23 of Ets results in a 14-fold decrease in k _cat/K _m, with little apparent change in k _cat. A peptide that binds to the DRS of ERK2 affects K _m, but not k _cat. Our kinetic analysis suggests that the unstructured N-terminus provides 10-fold uniform stabilization of the ground state ERK2•Ets•MgATP complex and intermediates of the enzymatic reaction.     

A bioreaction–diffusion model for growth of marine sponge explants in bioreactors

Marine sponges are sources of high-value bioactives. Engineering aspects of in vitro culture of sponges from cuttings (explants) are poorly understood. This work develops a diffusion-controlled growth model for sponge explants. The model assumes that the explant growth is controlled by diffusive transport of at least some nutrients from the surrounding medium into the explant that generally has a poorly developed aquiferous system for internal irrigation during early stages of growth. Growth is assumed to obey Monod-type kinetics. The model is shown to satisfactorily explain the measured growth behavior of the marine sponge Crambe crambe in two different growth media. In addition, the model is generally consistent with published data for growth of explants of the sponges Disidea avara and Hemimycale columella. The model predicted that nutrient concentration profiles for nutrients, such as dissolved oxygen within the explant, are consistent with data published by independent researchers. In view of the proposed model’s ability to explain available data for growth of several species of sponge explants, diffusive transport does play a controlling role in explant growth at least until a fully developed aquiferous system has become established. According to the model and experimental observations, the instantaneous growth rate depends on the size of the explant and all those factors that influence the diffusion of critical nutrients within the explant. Growth follows a hyperbolic profile that is consistent with the Monod kinetics.     

Hybrid MD-Nernst Planck model of α-hemolysin conductance properties

Motivated by experiments in which an applied electric field translocates polynucleotides through an α-hemolysin protein channel causing ionic current transient blockade, a hybrid simulation model is proposed to predict the conductance properties of the open channel. Time scales corresponding to ion permeation processes are reached using the Poisson–Nernst–Planck (PNP) electro-diffusion model in which both solvent and local ion concentrations are represented as a continuum. The diffusion coefficients of the ions (K⁺ and Cl^?) input in the PNP model are, however, calculated from all-atom molecular dynamics (MD). In the MD simulations, a reduced representation of the channel is used. The channel is solvated in a 1?M KCl solution, and an external electric field is applied. The pore specific diffusion coefficients for both ionic species are reduced 5–7 times in comparison to bulk values. Significant statistical variations (17–45%) of the pore-ions diffusivities are observed. Within the statistics, the ionic diffusivities remain invariable for a range of external applied voltages between 30 and 240?mV. In the 2D-PNP calculations, the pore stem is approximated by a smooth cylinder of radius ～9?Å with two constriction blocks where the radius is reduced to ～6?Å. The electrostatic potential includes the contribution from the atomistic charges. The MD-PNP model shows that the atomic charges are responsible for the rectifying behaviour and for the slight anion selectivity of the α-hemolysin pore. Independent of the hierarchy between the anion and cation diffusivities, the anionic contribution to the total ionic current will dominate. The predictions of the MD-PNP model are in good agreement with experimental data and give confidence in the present approach of bridging time scales by combining a microscopic and macroscopic model.     

A reaction–diffusion malaria model with incubation period in the vector population

Malaria is one of the most important parasitic infections in humans and more than two billion people are at risk every year. To understand how the spatial heterogeneity and extrinsic incubation period (EIP) of the parasite within the mosquito affect the dynamics of malaria epidemiology, we propose a nonlocal and time-delayed reaction–diffusion model. We then define the basic reproduction ratio R₀{\mathcal{R}_0} and show that R₀{\mathcal{R}_0} serves as a threshold parameter that predicts whether malaria will spread. Furthermore, a sufficient condition is obtained to guarantee that the disease will stabilize at a positive steady state eventually in the case where all the parameters are spatially independent. Numerically, we show that the use of the spatially averaged system may highly underestimate the malaria risk. The spatially heterogeneous framework in this paper can be used to design the spatial allocation of control resources.     

Spatial distribution of microalgae in marine systems: A reaction–diffusion model

In this paper, we have proposed a reaction–diffusion system of partial differential equations which model the plankton-nutrient interaction mediated by a toxin-determined functional response. It has been established that microalgae, a clean and green source of energy, can be potentially used for carbon capture and sequestration. The common biofuels (bio-diesel and ethanol) are efficiently extracted from microalgae of different shapes and sizes. A spatio-temporal model has been presented to guide exploration and harvesting of microalgae (e.g., dinoflagellates, cilliates, chlorella, etc.). The spatial distribution of the phytoplankton (microalgae) is determined by growth pattern of the biotic subsystem (phytoplankton and zooplankton); e.g., whether it is oscillatory or aperiodic. The model incorporates a toxin-determined functional response of the zooplankton, which can be parametrized for specific phytoplankton–zooplankton combinations in different aquatic bodies such as ponds, seas, and oceans. The present model does not take into account higher zooplankton’s role in maintaining the core subsystem. The temporal model is analytically investigated in terms of the existence criteria and stability analysis (both linear and nonlinear) of the possible equilibria and the spatio-temporal model is studied in terms of global stability, Turing instability and existence of Hopf-bifurcation which help us to explore the dynamical behavior of the spatial model system. Numerical simulations are carried out to support the obtained theoretical results. Simulation experiments and computed densities thereof (equal densities are codes by same color) suggest that the spatial distribution of microalgae is complex; e.g., spatial density of microalgae varies chaotically for certain parameter sets. Harvesting schedule can be designed based on information thus derived. It should be implemented carefully in case the spatial density distribution is chaotic.The sustainability of the marine system for future use has been the prime concern. Parameters of harvesting strategy (time, intensity and technology) are determined in such a way that exploitation causes minimal damage to the environment and the yield of the harvest is maximal. Future studies would consider larger carnivorous fishes (e.g., Squids, Dolphins) on system’s dynamics. The effect of oceanic noise and colloidal swarming of zooplankton in the presence of bacteria will also be incorporated.     

Origins and suppression of oscillations in a computational model of Parkinson’s disease

Efficacy of deep brain stimulation (DBS) for motor signs of Parkinson’s disease (PD) depends in part on post-operative programming of stimulus parameters. There is a need for a systematic approach to tuning parameters based on patient physiology. We used a physiologically realistic computational model of the basal ganglia network to investigate the emergence of a 34 Hz oscillation in the PD state and its optimal suppression with DBS. Discrete time transfer functions were fit to post-stimulus time histograms (PSTHs) collected in open-loop, by simulating the pharmacological block of synaptic connections, to describe the behavior of the basal ganglia nuclei. These functions were then connected to create a mean-field model of the closed-loop system, which was analyzed to determine the origin of the emergent 34 Hz pathological oscillation. This analysis determined that the oscillation could emerge from the coupling between the globus pallidus external (GPe) and subthalamic nucleus (STN). When coupled, the two resonate with each other in the PD state but not in the healthy state. By characterizing how this oscillation is affected by subthreshold DBS pulses, we hypothesize that it is possible to predict stimulus frequencies capable of suppressing this oscillation. To characterize the response to the stimulus, we developed a new method for estimating phase response curves (PRCs) from population data. Using the population PRC we were able to predict frequencies that enhance and suppress the 34 Hz pathological oscillation. This provides a systematic approach to tuning DBS frequencies and could enable closed-loop tuning of stimulation parameters.     

Instability in a generalized Keller–Segel model

We present a generalized Keller–Segel model where an arbitrary number of chemical compounds react, some of which are produced by a species, and one of which is a chemoattractant for the species. To investigate the stability of homogeneous stationary states of this generalized model, we consider the eigenvalues of a linearized system. We are able to reduce this infinite dimensional eigenproblem to a parametrized finite dimensional eigenproblem. By matrix theoretic tools, we then provide easily verifiable sufficient conditions for destabilizing the homogeneous stationary states. In particular, one of the sufficient conditions is that the chemotactic feedback is sufficiently strong. Although this mechanism was already known to exist in the original Keller–Segel model, here we show that it is more generally applicable by significantly enlarging the class of models exhibiting this instability phenomenon which may lead to pattern formation.