The motion of a eukaryotic cell presents a variety of interesting and challenging problems from both a modeling and a computational perspective. The processes span many spatial scales (from molecular to tissue) as well as disparate time scales, with reaction kinetics on the order of seconds, and the deformation and motion of the cell occurring on the order of minutes. The computational difficulty, even in 2D, resides in the fact that the problem is inherently one of deforming, non-stationary domains, bounded by an elastic perimeter, inside of which there is redistribution of biochemical signaling substances. Here we report the results of a computational scheme using the immersed boundary method to address this problem. We adopt a simple reaction-diffusion system that represents an internal regulatory mechanism controlling the polarization of a cell, and determining the strength of protrusion forces at the front of its elastic perimeter. Using this computational scheme we are able to study the effect of protrusive and elastic forces on cell shapes on their own, the distribution of the reaction-diffusion system in irregular domains on its own, and the coupled mechanical-chemical system. We find that this representation of cell crawling can recover important aspects of the spontaneous polarization and motion of certain types of crawling cells.

A numerical model based on one-dimensional balance laws and ad hoc zero-dimensional boundary conditions is tested against experimental data. The study concentrates on the circle of Willis, a vital subnetwork of the cerebral vasculature. The main goal is to obtain efficient and reliable numerical tools with predictive capabilities. The flow is assumed to obey the Navier-Stokes equations, while the mechanical reactions of the arterial walls follow a viscoelastic model. Like many previous studies, a dimension reduction is performed through averaging. Unlike most previous work, the resulting model is both calibrated and validated against in vivo data, more precisely transcranial Doppler data of cerebral blood velocity. The network considered has three inflow vessels and six outflow vessels. Inflow conditions come from the data, while outflow conditions are modeled. Parameters in the outflow conditions are calibrated using a subset of the data through ensemble Kalman filtering techniques. The rest of the data is used for validation. The results demonstrate the viability of the proposed approach.

Time series associated with single-molecule experiments and/or simulations contain a wealth of multiscale information about complex biomolecular systems. We demonstrate how a collection of Penalized-splines (P-splines) can be useful in quantitatively summarizing such data. In this work, functions estimated using P-splines are associated with stochastic differential equations (SDEs). It is shown how quantities estimated in a single SDE summarize fast-scale phenomena, whereas variation between curves associated with different SDEs partially reflects noise induced by motion evolving on a slower time scale. P-splines assist in "semiparametrically" estimating nonlinear SDEs in situations where a time-dependent external force is applied to a single-molecule system. The P-splines introduced simultaneously use function and derivative scatterplot information to refine curve estimates. We refer to the approach as the PuDI (P-splines using Derivative Information) method. It is shown how generalized least squares ideas fit seamlessly into the PuDI method. Applications demonstrating how utilizing uncertainty information/approximations along with generalized least squares techniques improve PuDI fits are presented. Although the primary application here is in estimating nonlinear SDEs, the PuDI method is applicable to situations where both unbiased function and derivative estimates are available.

We present a new (DENM) that describes the unfolding process of a force-loaded protein. The protein interaction network and its potentials are constructed based on information of its native-state structure obtained from the Protein Data Bank, with network nodes positioned at the coordinates of the protein backbone. Specifically, to mimic the unfolding process, i.e., to simulate the process of overcoming the local energy barrier on the free energy landscape with force loading, the noncovalent protein network bonds (i.e., hydrogen bonds, salt bridges, hydrophobic contacts, etc.) are broken one-by-one with a certain probability, while the strong covalent bonds along the backbone (i.e., peptide bonds, disulfide bonds, etc.) are kept intact. The jumping event from local energy minima (bonds breaking rate) are chosen according to Kramer's theory and the Bell model. Moreover, we exploit the self-similar structure of proteins at different scales to design an effective coarse-graining procedure for DENM with optimal parameter selection. The robustness of DENM is validated by coarse-grained molecular dynamics (MD) simulation against atomistic MD simulation of force-extension processes of the Fibrinogen and Titin Immunoglobulin proteins. We observe that the native structure of the proteins determines the unfolding dynamics (including large deviations) and not just the fluctuations around the native state.

Subdiffusion has been proposed as an explanation of various kinetic phenomena inside living cells. In order to fascilitate large-scale computational studies of subdiffusive chemical processes, we extend a recently suggested mesoscopic model of subdiffusion into an accurate and consistent reaction-subdiffusion computational framework. Two different possible models of chemical reaction are revealed and some basic dynamic properties are derived. In certain cases those mesoscopic models have a direct interpretation at the macroscopic level as fractional partial differential equations in a bounded time interval. Through analysis and numerical experiments we estimate the macroscopic effects of reactions under subdiffusive mixing. The models display properties observed also in experiments: for a short time interval the behavior of the diffusion and the reaction is ordinary, in an intermediate interval the behavior is anomalous, and at long times the behavior is ordinary again.

Biomolecules exhibit conformational fluctuations near equilibrium states, inducing uncertainty in various biological properties in a way. We have developed a general method to quantify the uncertainty of target properties induced by conformational fluctuations. Using a generalized polynomial chaos (gPC) expansion, we construct a surrogate model of the target property with respect to varying conformational states. To alleviate the high-dimensionality of the corresponding stochastic space, we propose a method to increase the sparsity of the gPC expansion by defining a set of conformational "active space" random variables. With the increased sparsity, we employ the compressive sensing method to accurately construct the surrogate model. We demonstrate the performance of the surrogate model by evaluating fluctuation-induced uncertainty in solvent-accessible surface area for the bovine trypsin inhibitor protein system and show that the new approach offers more accurate statistical information than standard Monte Carlo approaches. Furthermore, the constructed surrogate model also enables us to evaluate the target property under various conformational states, yielding a more accurate response surface than standard sparse grid collocation methods. In particular, the new method provides higher accuracy in high-dimensional systems, such as biomolecules, where sparse grid performance is limited by the accuracy of the computed quantity of interest. Our new framework is generalizable and can be used to investigate the uncertainty of a wide variety of target properties in biomolecular systems.

The distinct timescales of synaptic plasticity and neural activity dynamics play an important role in the brain's learning and memory systems. Activity-dependent plasticity reshapes neural circuit architecture, determining spontaneous and stimulus-encoding spatiotemporal patterns of neural activity. Neural activity bumps maintain short term memories of continuous parameter values, emerging in spatially organized models with short-range excitation and long-range inhibition. Previously, we demonstrated nonlinear Langevin equations derived using an interface method which accurately describe the dynamics of bumps in continuum neural fields with separate excitatory/inhibitory populations. Here we extend this analysis to incorporate effects of short term plasticity that dynamically modifies connectivity described by an integral kernel. Linear stability analysis adapted to these piecewise smooth models with Heaviside firing rates further indicates how plasticity shapes the bumps' local dynamics. Facilitation (depression), which strengthens (weakens) synaptic connectivity originating from active neurons, tends to increase (decrease) stability of bumps when acting on excitatory synapses. The relationship is inverted when plasticity acts on inhibitory synapses. Multiscale approximations of the stochastic dynamics of bumps perturbed by weak noise reveal that the plasticity variables evolve to slowly diffusing and versions of their stationary profiles. Nonlinear Langevin equations associated with bump positions or interfaces coupled to slowly evolving projections of plasticity variables accurately describe how these smoothed synaptic efficacy profiles can tether or repel wandering bumps.

A newly developed, coarse-grained treatment of the low-frequency normal modes of DNA has been adapted to study the torsional properties of fully extended, double-helical molecules. Each base pair is approximated in this scheme as a rigid body, and molecular structure is described in terms of the relative position and orientation of successive base pairs. The torsional modulus C is computed from the lowest-frequency normal twisting mode using expressions valid for a homogeneous, naturally straight elastic rod. Fluctuations of local dimeric structure, including the coupled variation of conformational parameters, are based on the observed arrangements of neighboring base pairs in high-resolution structures. Chain ends are restrained by an elastic energy term. The calculations show how the end-to-end constraints placed on a naturally straight DNA molecule, in combination with the natural conformational features of the double helix, can account for the substantially larger torsional moduli determined with state-of-the-art, single-molecule experiments compared to values extracted from solution measurements and/or incorporated into theories to account for the force-extension properties of single molecules. The computed normal-mode frequencies and torsional moduli increase substantially if base pairs are inclined with respect to the double-helical axis and the deformations of selected conformational variables follow known interdependent patterns. The changes are greatest if the fluctuations in dimeric twisting are coupled with parameters that directly alter the end-to-end displacement. Imposed restraints that mimic the end-to-end conditions of single-molecule experiments then impede the twisting of base pairs and increase the torsional modulus. The natural inclination of base pairs concomitantly softens the Young's modulus, i.e., ease of duplex stretching. The analysis of naturally curved DNA points to a drop in the torsional modulus upon imposed extension of the double-helical molecule.

Numerous centrality measures have been developed to quantify the importances of nodes in time-independent networks, and many of them can be expressed as the leading eigenvector of some matrix. With the increasing availability of network data that changes in time, it is important to extend such eigenvector-based centrality measures to time-dependent networks. In this paper, we introduce a principled generalization of network centrality measures that is valid for any eigenvector-based centrality. We consider a temporal network with nodes as a sequence of layers that describe the network during different time windows, and we couple centrality matrices for the layers into a matrix of size × whose dominant eigenvector gives the centrality of each node at each time . We refer to this eigenvector and its components as a , as it reflects the importances of both the node and the time layer . We also introduce the concepts of and centralities, which facilitate the study of centrality trajectories over time. We find that the strength of coupling between layers is important for determining multiscale properties of centrality, such as localization phenomena and the time scale of centrality changes. In the strong-coupling regime, we derive expressions for , which are given by the zeroth-order terms of a singular perturbation expansion. We also study first-order terms to obtain , which concisely describe the magnitude of nodes' centrality changes over time. As examples, we apply our method to three empirical temporal networks: the United States Ph.D. exchange in mathematics, costarring relationships among top-billed actors during the Golden Age of Hollywood, and citations of decisions from the United States Supreme Court.

Using particle-scale models to accurately describe property enhancements and phase transitions in macroscopic behavior is a major engineering challenge in composite materials science. To address some of these challenges, we use the graph theoretic property of rigidity to model mechanical reinforcement in composites with stiff rod-like particles. We develop an efficient algorithmic approach called (RGC) to describe the transition from floppy to rigid in disordered fiber networks ("rod-hinge systems"), which form the reinforcing phase in many composite systems. To establish RGC on a firm theoretical foundation, we adapt rigidity matroid theory to identify primitive topological network motifs that serve as rules for composing interacting rigid particles into larger rigid components. This approach is computationally efficient and stable, because RGC requires only topological information about rod interactions (encoded by a sparse unweighted network) rather than geometrical details such as rod locations or pairwise distances (as required in rigidity matroid theory). We conduct numerical experiments on simulated two-dimensional rod-hinge systems to demonstrate that RGC closely approximates the rigidity percolation threshold for such systems, through comparison with the pebble game algorithm (which is exact in two dimensions). Importantly, whereas the pebble game is derived from Laman's condition and is only valid in two dimensions, the RGC approach naturally extends to higher dimensions.

We present the first mathematical model of flow-mediated primary hemostasis in an extravascular injury which can track the process from initial deposition to occlusion. The model consists of a system of ordinary differential equations (ODEs) that describe platelet aggregation (adhesion and cohesion), soluble-agonist-dependent platelet activation, and the flow of blood through the injury. The formation of platelet aggregates increases resistance to flow through the injury, which is modeled using the Stokes-Brinkman equations. Data from analogous experimental (microfluidic flow) and partial differential equation models informed parameter values used in the ODE model description of platelet adhesion, cohesion, and activation. This model predicts injury occlusion under a range of flow and platelet activation conditions. Simulations testing the effects of shear and activation rates resulted in delayed occlusion and aggregate heterogeneity. These results validate our hypothesis that flow-mediated dilution of activating chemical adenosine diphosphate hinders aggregate development. This novel modeling framework can be extended to include more mechanisms of platelet activation as well as the addition of the biochemical reactions of coagulation, resulting in a computationally efficient high throughput screening tool of primary and secondary hemostasis.

We propose parameter optimization techniques for weighted ensemble sampling of Markov chains in the steady-state regime. Weighted ensemble consists of replicas of a Markov chain, each carrying a weight, that are periodically resampled according to their weights inside of each of a number of bins that partition state space. We derive, from first principles, strategies for optimizing the choices of weighted ensemble parameters, in particular the choice of bins and the number of replicas to maintain in each bin. In a simple numerical example, we compare our new strategies with more traditional ones and with direct Monte Carlo.

The mechanism of gas transport across cell membranes remains a topic of considerable interest, particularly regarding the extent to which lipids vs. specific membrane proteins provide conduction pathways. Studies of transmembrane (CO) transport often rely on data collected under controlled conditions, using pH-sensitive microelectrodes at the extracellular surface to record changes due to extracellular CO diffusion and reactions. Although recent detailed computational models can predict a qualitatively correct behavior, a mismatch between the dynamical ranges of the predicted and observed pH curves raises the question whether the discrepancy may be due to a bias introduced by the pH electrode itself. More specifically, it is reasonable to ask whether bringing the electrode tip near or in contact with the membrane creates a local microenvironment between the electrode tip and the membrane, so that the measured data refer to the microenvironment rather than to the free surface. Here, we introduce a detailed computational model, designed to address this question. We find that, as long as a zone of free diffusion exists between the tip and the membrane, the microenvironment behaves effectively as the free membrane. However, according to our model, when the tip contacts the membrane, partial quenching of extracellular diffusion by the electrode rim leads to a significant increase in the pH dynamics under the electrode, matching values measured in physiological experiments. The computational schemes for the model predictions are based on semi-discretization by a finite-element method, and an implicit-explicit time integration scheme to capture the different time scales of the system.

We present a novel weak formulation and discretization for discovering governing equations from noisy measurement data. This method of learning differential equations from data fits into a new class of algorithms that replace pointwise derivative approximations with linear transformations and variance reduction techniques. Compared to the standard SINDy algorithm presented in [S. L. Brunton, J. L. Proctor, and J. N. Kutz, , 113 (2016), pp. 3932-3937], our so-called weak SINDy (WSINDy) algorithm allows for reliable model identification from data with large noise (often with ratios greater than 0.1) and reduces the error in the recovered coefficients to enable accurate prediction. Moreover, the coefficient error scales linearly with the noise level, leading to high-accuracy recovery in the low-noise regime. Altogether, WSINDy combines the simplicity and efficiency of the SINDy algorithm with the natural noise reduction of integration, as demonstrated in [H. Schaeffer and S. G. McCalla, , 96 (2017), 023302], to arrive at a robust and accurate method of sparse recovery.

The discrete chemical master equation (dCME) provides a fundamental framework for studying stochasticity in mesoscopic networks. Because of the multi-scale nature of many networks where reaction rates have large disparity, directly solving dCMEs is intractable due to the exploding size of the state space. It is important to truncate the state space effectively with quantified errors, so accurate solutions can be computed. It is also important to know if all major probabilistic peaks have been computed. Here we introduce the Accurate CME (ACME) algorithm for obtaining direct solutions to dCMEs. With multi-finite buffers for reducing the state space by (!), exact steady-state and time-evolving network probability landscapes can be computed. We further describe a theoretical framework of aggregating microstates into a smaller number of macrostates by decomposing a network into independent aggregated birth and death processes, and give an method for rapidly determining steady-state truncation errors. The maximal sizes of the finite buffers for a given error tolerance can also be pre-computed without costly trial solutions of dCMEs. We show exactly computed probability landscapes of three multi-scale networks, namely, a 6-node toggle switch, 11-node phage-lambda epigenetic circuit, and 16-node MAPK cascade network, the latter two with no known solutions. We also show how probabilities of rare events can be computed from first-passage times, another class of unsolved problems challenging for simulation-based techniques due to large separations in time scales. Overall, the ACME method enables accurate and efficient solutions of the dCME for a large class of networks.

In this paper we present a method for the multi-resolution comparison of biomolecular electrostatic potentials without the need for global structural alignment of the biomolecules. The underlying computational geometry algorithm uses multi-resolution attributed contour trees (MACTs) to compare the topological features of volumetric scalar fields. We apply the MACTs to compute electrostatic similarity metrics for a large set of protein chains with varying degrees of sequence, structure, and function similarity. For calibration, we also compute similarity metrics for these chains by a more traditional approach based upon 3D structural alignment and analysis of Carbo similarity indices. Moreover, because the MACT approach does not rely upon pairwise structural alignment, its accuracy and efficiency promises to perform well on future large-scale classification efforts across groups of structurally-diverse proteins. The MACT method discriminates between protein chains at a level comparable to the Carbo similarity index method; i.e., it is able to accurately cluster proteins into functionally-relevant groups which demonstrate strong dependence on ligand binding sites. The results of the analyses are available from the linked web databases http://ccvweb.cres.utexas.edu/MolSignature/ and http://agave.wustl.edu/similarity/. The MACT analysis tools are available as part of the public domain library of the Topological Analysis and Quantitative Tools (TAQT) from the Center of Computational Visualization, at the University of Texas at Austin (http://ccvweb.csres.utexas.edu/software). The Carbo software is available for download with the open-source APBS software package at http://apbs.sf.net/.

We give a mathematical framework for Exact Milestoning, a recently introduced algorithm for mapping a continuous time stochastic process into a Markov chain or semi-Markov process that can be efficiently simulated and analyzed. We generalize the setting of Exact Milestoning and give explicit error bounds for the error in the Milestoning equation for mean first passage times.

Cell crawling is an important biological phenomenon underlying coordinated cell movement in morphogenesis, cancer, and wound healing. In recent decades the process of cell crawling has been experimentally and theoretically dissected into further subprocesses: protrusion of the cell at its leading edge, retraction of the cell body, and graded adhesion. A number of one-dimensional (1-D) models explain successfully a proximal-distal organization and movement of the motile cell. However, more adequate two-dimensional (2-D) models are lacking. We propose a multiscale 2-D computational model of the lamellipodium (motile appendage) of a simply shaped, rapidly crawling fish keratocyte cell. We couple submodels of (i) protrusion and adhesion at the leading edge, (ii) the elastic 2-D lamellipodial actin network, (iii) the actin-myosin contractile bundle at the rear edge, and (iv) the convection-reaction-diffusion actin transport on the free boundary lamellipodial domain. We simulate the combined model numerically using a finite element approach. The simulations reproduce observed cell shapes, forces, and movements and explain some experimental results on perturbations of the actin machinery. This novel 2-D model of the crawling cell makes testable predictions and posits questions to be answered by future modeling.

In biomechanics, local phenomena, such as tissue perfusion, are strictly related to the global features of the surrounding blood circulation. In this paper, we propose a heterogeneous model where a local, accurate, 3D description of tissue perfusion by means of fluid flows through deformable porous media equations is coupled with a systemic, 0D, lumped model of the remainder of the circulation, where the fluid flow through a vascular network is described via its analog with a current flowing through an electric circuit. This represents a multiscale strategy, which couples an initial boundary value problem to be used in a specific tissue region with an initial value problem in the surrounding circulatory system. This PDE/ODE coupling leads to interface conditions enforcing the continuity of mass and the balance of stresses across models at different scales, and careful consideration is taken to address this interface mismatch. The resulting system involves PDEs of mixed type with interface conditions depending on nonlinear ODEs. A new result on local existence of solutions for this multiscale interface coupling is provided in this article.