Past OCCAM Wednesday morning Events

Russian mathematician V.I.Arnold conjectured that convex, homogeneous bodies with less than four equilibria (also called mono-monostatic) may exist. Not only did his conjecture turn out to be true, the newly discovered objects show various interesting features. Our goal is to give an overview of these findings as well as to present some new results. We will point out that mono-monostatic bodies are neither flat, nor thin, they are not similar to typical objects with more equilibria and they are hard to approximate by polyhedra. Despite these "negative" traits, there seems to be strong indication that these forms appear in Nature due to their special mechanical properties.

Telomeres, non-coding terminal structures of DNA strands, consist of repetitive long tandem repeats of a specific length. An absence of an enzyme, telomerase, in certain cellular structures requires an alternative telomerase-independent pathway for telomeric sequence length regulation. Besides linear telomeres other configurations such as telomeric circles and telomeric loops were experimentally observed. They are suspected to play an important role in a universal mechanism for stabilization of the ends of linear DNA that possibly dates back to pre-telomerase ages. We propose a mathematical model that captures biophysical interactions of various telomeric structures on a short time scale and that is able to reproduce experimental measurements in mtDNA of yeast. Moreover, the model opens up a couple of interesting mathematical problems such as validity of a quasi-steady state approximation and dynamic properties of discrete coagulation-fragmentation systems. We also identify and estimate key factors influencing the length distribution of telomeric circles, loops and strand invasions using numerical simulations.

For a given positive measure on a fixed domain, Gaussian quadrature routines can be defined via their property of integrating an arbitrary polynomial of degree $2n+1$ exactly using only $n+1$ quadrature nodes. In the special case of Gauss--Jacobi quadrature, this means that $$\int_{-1}^1 (1+x)^\alpha(1-x)^\beta f(x) dx = \sum_{j=0}^{n} w_j f(x_j), \quad \alpha, \beta > -1, $$ whenever $f(x)$ is a polynomial of degree at most $2n+1$. When $f$ is not a polynomial, but a function analytic in a neighbourhood of $[-1,1]$, the above is not an equality but an approximation that converges exponentially fast as $n$ is increased.

An undergraduate mathematician taking a numerical analysis course could tell you that the nodes $x_j$ are roots of the Jacobi polynomial $P^{\alpha,\beta}_{n+1}(x)$, the degree $n+1$ polynomial orthogonal with respect to the given weight, and that the quadrature weight at each node is related to the derivative $P'^{\alpha,\beta}_{n+1}(x_j)$. However, these values are not generally known in closed form, and we must compute them numerically... but how?

Traditional approaches involve applying the recurrence relation satisfied by the orthogonal polynomials, or solving the Jacobi matrix eigenvalue problem in the algorithm of Golub and Welsch, but these methods are inherently limited by a minimal complexity of $O(n^2)$. The current state-of-the-art is the $O(n)$ algorithm developed by Glasier, Liu, and Rokhlin, which hops from root to root using a Newton iteration evaluated with a Taylor series defined by the ODE satisfied by $P^{\alpha,\beta}_{n+1}$.

We propose an alternative approach, whereby the function and derivative evaluations required in the Newton iteration are computed independently in $O(1)$ operations per point using certain well-known asymptotic expansions. We shall review these expansions, outline the new algorithm, and demonstrate improvements in both accuracy and efficiency. 

We present some recent results concerning domain wall motion in one-dimensional nanowires, including the existence, velocity and stability of travelling-wave and precessing solutions.  We consider the case of unixial anisotropy, characteristic of wires with symmetrical (e.g., circular) cross-section, as opposed to strongly anisotropic geometries (films and strips) that have received greater attention.  This is joint work with Arseni Goussev and Valeriy Slastikov.

Actin polymerization in vivo is regulated spatially and temporally by a web of signalling proteins. We developed detailed physico-chemical, stochastic models of lamellipodia and filopodia, which are projected by eukaryotic cells during cell migration, and contain dynamically remodelling actin meshes and bundles. In a recent work we studied how molecular motors regulate growth dynamics of elongated organelles of living cells. We determined spatial distributions of motors in such organelles, corresponding to a basic scenario when motors only walk along the substrate, bind, unbind, and diffuse. We developed a mean field model, which quantitatively reproduces elaborate stochastic simulation results as well as provides a physical interpretation of experimentally observed distributions of Myosin IIIa in stereocilia and filopodia. The mean field model showed that the jamming of the walking motors is conspicuous, and therefore damps the active motor flux. However, when the motor distributions are coupled to the delivery of actin monomers towards the tip, even the concentration bump of G-actin that they create before they jam is enough to speed up the diffusion to allow for severalfold longer filopodia. We found that the concentration profile of G-actin along the filopodium is rather non-trivial, containing a narrow minimum near the base followed by a broad maximum. For efficient enough actin transport, this non-monotonous shape is expected to occur under a broad set of conditions. We also find that the stationary motor distribution is universal for the given set of model parameters regardless of the organelle length, which follows from the form of the kinetic equations and the boundary conditions.

Particle-based stochastic reaction-diffusion models have recently been used to study a number of problems in cell biology. These methods are of interest when both noise in the chemical reaction process and the explicit motion of molecules are important. Several different mathematical models have been used, some spatially-continuous and others lattice-based. In the former molecules usually move by Brownian Motion, and may react when approaching each other. For the latter molecules undergo continuous time random-walks, and usually react with fixed probabilities per unit time when located at the same lattice site.

As motivation, we will begin with a brief discussion of the types of biological problems we are studying and how we have used stochastic reaction-diffusion models to gain insight into these systems. We will then introduce several of the stochastic reaction-diffusion models, including the spatially continuous Smoluchowski diffusion limited reaction model and the lattice-based reaction-diffusion master equation. Our work studying the rigorous relationships between these models will be presented. Time permitting, we may also discuss some of our efforts to develop improved numerical methods for solving several of the models.

The reaction-diffusion master equation (RDME) is a popular model in systems biology. In the RDME, diffusion is modeled as discrete jumps between voxels in the computational domain. However, it has been demonstrated that a more fine-grained model is required to resolve all the dynamics of some highly diffusion-limited systems.

In Greenʼs Function Reaction Dynamics (GFRD), a method based on the Smoluchowski model, diffusion is modeled continuously in space.

This will be more accurate at fine scales, but also less efficient than a discrete-space model.

We have developed a hybrid method, combining the RDME and the GFRD method, making it possible to do the more expensive fine-grained simulations only for the species and in the parts of space where it is required in order to resolve all the dynamics, and more coarse-grained simulations where possible. We have applied this method to a model of a MAPK-pathway, and managed to reduce the number of molecules simulated with GFRD by orders of magnitude and without an appreciable loss of accuracy.

In this talk, I will introduce stochastic models to describe the state of the chemical networks using continuous-time Markov chains.
First, I will talk about the multiscale approximation method developed by Ball, Kurtz, Popovic, and Rempala (2006). Extending their method, we construct a general multiscale approximation in chemical reaction networks. We embed a stochastic model for a chemical reaction network into a family of models parameterized by a large parameter N. If reaction rate constants and species numbers vary over a wide range, we scale these numbers by powers of the parameter N. We develop a systematic approach to choose an appropriate set of scaling exponents. When the scaling suggests subnetworks have di erent time-scales, the subnetwork in each time scale is approximated by a limiting model involving a subset of reactions and species.

After that, I will briefly introduce Gaussian approximation using a central limit theorem, which gives a model with more detailed uctuations which may be not captured by the limiting models in multiscale approximations.

Next, we consider modeling of a chemical network with both reaction and diffusion.
We discretize the spatial domain into several computational cells and model diffusion as a reaction where the molecule of species in one computational cell moves to the neighboring one. In this case, the important question is how many computational cells we need to use for discretization to get balance between e ective diffusion rates and reaction rates both of which depend on the computational cell size. We derive a condition under which concentration of species converges to its uniform solution exponentially. Replacing a system domain size in this condition by computational cell size in our stochastic model, we derive an upper bound
for the computational cell size.

Finally, I will talk about stochastic reaction-diffusion models of pattern formation. Spatially distributed signals called morphogens influence gene expression that determines phenotype identity of cells. Generally, different cell types are segregated by boundary between
them determined by a threshold value of some state variables. Our question is how sensitive the location of the boundary to variation in parameters. We suggest a stochastic model for boundary determination using signaling schemes for patterning and investigate their effects on the variability of the boundary determination between cells.

Atherosclerosis is the leading cause of death, both above and below age 65, in the United States and all Western countries. Its earliest prelesion events appear to be the transmural (across the wall)-pressure (DP)-driven advection of large molecules such as low-density lipoprotein (LDL) cholesterol from the blood into the inner wall layers across the monolayer of endothelial cells that tile the blood-wall interface. This transport occurs through the junctions around rare (~one cell every few thousand) endothelial cells whose junctions are wide enough to allow large molecules to pass. These LDL molecules can bind to extracellular matrix (ECM) in the wall’s thin subendothelial intima (SI) layer and accumulate there. On the other hand, the overall transmural water flow can dilute the local intima LDL concentration, thereby slowing its kinetics of binding to ECM, and flushes unbound lipid from the wall. An understanding of the nature of this water flow is clearly critical.

            We have found that rat aortic endothelial cells express the ubiquitous membrane water-channel protein aquaporin-1 (AQP), and that blocking its water channel or knocking down its expression significantly reduces the apparent hydraulic conductivity Lp of the endothelium and, consequently of the entire wall. This decrease has an unexpected and strong DP -dependence. We present a fluid mechanics theory based on the premise that DP compacts the SI, which, as we show, lowers its Lp. The theory shows that blocking or knocking down AQP flow changes the critical DP at which this compaction occurs and explains our observed dependence of Lp on DP. Such compaction may affect lipid transport and accumulation in vivo. However, AQP’s sharp water selectivity gives rise to an oncotic paradox: the SI should quickly become hypotonic and shut down this AQP flow. The mass transfer problem resolve this paradox. The importance of aquaporin-based, rather than simply junctional water transport is that transport via protein channels allows for the possibility of active control of vessel Lp by up- or down-regulation of protein expression. We show that rat aortic endothelial cells significantly change their AQP numbers in response to chronic hypertension (high blood pressure), which may help explain the as yet poorly-understood fact that hypertension correlates with atherosclerosis. We also consider lowering AQP numbers as a strategy to affect disease progression.

Kernel functions are suitable tools for multivariate scattered data approximation. In this talk, we discuss the conditioning and stability of optimal reconstruction schemes from multivariate scattered data by using

(conditionally) positive definite kernel functions. Our discussion first provides basic Riesz-type stability estimates for the utilized reconstruction method, before particular focus is placed on upper and lower bounds of the Lebesgue constants.

If time allows, we will finally draw our attention to relevant aspects concerning the stability of penalized least squares approximation.

In many fields of science and engineering, such as fluid or structural mechanics and electric circuit design, large scale dynamical systems need to be simulated, optimized or controlled. They are often described by discretizations of systems of nonlinear partial differential equations yielding high-dimensional discrete phase spaces. For this reason, in recent decades, research has mainly focused on the development of sophisticated analytical and numerical tools to help understand the overall system behavior. During this time meshless methods have enjoyed significant interest in the research community and in some commercial simulators (e.g., LS-DYNA). In this talk I will describe a normalized-corrected meshless method which ensures linear completeness and improved accuracy. The resulting scheme not only provides first order consistency O(h) but also alleviates the particle deficiency (kernel support incompleteness) problem at the boundary. Furthermore, a number of improvements to the kernel derivative approximation are proposed.

To illustrate the performance of the meshless method, we present results for different problems from various fields of science and engineering (i.e. nano-tubes modelling, solid mechanics, damage mechanics, fluid mechanics, coupled interactions of solids and fluids). Special attention is paid to fluid flow in porous media. The primary attraction of the present approach is that it provides a weak formulation for Darcy's law which can be used in further development of meshless methods.

The graph realization problem has received a great deal of attention in recent years, due to its importance in applications such as wireless sensor networks and structural biology. We introduce the ASAP algorithm, for the graph realization problem in R^d, given a sparse and noisy set of distance measurements associated to the edges of a globally rigid graph. ASAP is a divide and conquer, non-incremental and non-iterative algorithm, which integrates local distance information into a global structure determination. Our approach starts with identifying, for every node, a subgraph of its 1-hop neighborhood graph, which can be accurately embedded in its own coordinate system. In the noise-free case, the computed coordinates of the sensors in each patch must agree with their global positioning up to some unknown rigid motion, that is, up to translation, rotation and possibly reflection. In other words, to every patch there corresponds an element of the Euclidean group Euc(3) of rigid transformations in R^3, and the goal is to estimate the group elements that will properly align all the patches in a globally consistent way. The reflections and rotations are estimated using a recently developed eigenvector synchronization algorithm, while the translations are estimated by solving an overdetermined linear system. Furthermore, the algorithm successfully incorporates information specific to the molecule problem in structural biology, in particular information on known substructures and their orientation. In addition, we also propose SP-ASAP, a faster version of ASAP, which uses a spectral partitioning algorithm as a preprocessing step for dividing the initial graph into smaller subgraphs. Our extensive numerical simulations show that ASAP and SP-ASAP are very robust to high levels of noise in the measured distances and to sparse connectivity in the measurement graph, and compare favorably to similar state-of-the art localization algorithms. Time permitting, we briefly discuss the analogy between the graph realization and the low-rank matrix completion problems, as well as an application of synchronization over Z_2 and its variations to bipartite multislice networks.

Mathematical models of chemotactic movement of bacterial populations are often written as systems of partial differential equations for the densities of bacteria and concentrations of extracellular signaling chemicals. This approach has been employed since the seminal work of Keller and Segel in the 1970s [Keller and Segel, J. Theor. Biol., 1971]. The system has been shown to permit travelling wave solutions which correspond to travelling band formation in bacterial colonies, yet only under specific criteria, such as a singularity in the chemotactic sensitivity function as the signal approaches zero. Such a singularity generates infinite macroscopic velocities that ar biologically unrealistic. Here we present a microscopic model that takes into consideration relevant details of the intracellular processes while avoiding the singularity in the chemotactic sensitivity. We show that this model permits travelling wave solutions and predicts the formation of other bacterial patterns such as radial and spiral streams. We also present connections of this microscopic model with macroscopic models of bacterial chemotaxis. This is joint work with Radek Erban, Benjamin Franz, Hyung Ju Hwang, and Kevin J.

Painter.

Many swimming microorganisms inhabit heterogeneous environments consisting of solid particles immersed in viscous fluid. Such environments require the organisms attempting to move through them to negotiate both hydrodynamic forces and geometric constraints. Here, we study this kind of locomotion by first observing the kinematics of the small nematode and model organism Caenorhabditis elegans in fluid-filled, micro-pillar arrays. We then compare its dynamics with those given by numerical simulations of a purely mechanical worm model that accounts only for the hydrodynamic and contact interactions with the obstacles. We demonstrate that these interactions allow simple undulators to achieve speeds as much as an order of magnitude greater than their free-swimming values. More generally, what appears as behavior and sensing can sometimes be explained through simple mechanics.

We will consider a simplified model for on-chip power distribution networks of array bonded integrated circuits. In this model the voltage is the solution of a Poisson equation in an infinite planar domain whose boundary is an array of circular or square pads of size $\epsilon$. We deal with the singular limit as $\epsilon\to 0$ and we are interested in deriving an explicit formula for the maximum voltage drop in the domain in terms of a power series in $\epsilon$. A procedure based on the method of matched asymptotic expansions will be presented to compute all the successive terms in the approximation, which can be interpreted as using multipole solutions of equations involving spatial derivatives of $\delta$-functions.

Fractional differential equations are becoming increasingly used as a modelling tool for processes associated with anomalous diffusion or spatial heterogeneity. However, the presence of a fractional differential operator causes memory (time fractional) or nonlocality (space fractional) issues that impose a number of computational constraints. In this talk we discuss efficient, scalable techniques for solving fractional-in-space reaction diffusion equations combining the finite element method with robust techniques for computing the fractional power of a matrix times a vector. We shall demonstrate the methods on a number examples which show the qualitative difference in solution profiles between standard and fractional diffusion models.

Diffusive process with discontinuous coefficients provide significant computational challenges. We consider the solution of a diffusive process in a domain where the diffusion coefficient changes discontinuously across a curved interface. Rather than seeking to construct discretizations that match the interface, we consider the use of regularly-shaped meshes so that the interface "cuts'' through the cells (elements or volumes). Consequently, the discontinuity in the diffusion coefficients has a strong impact on the accuracy and convergence of the numerical method. We develop an adjoint based a posteriori error analysis technique to estimate the error in a given quantity of interest (functional of the solution). In order to employ this method, we first construct a systematic approach to discretizing a cut-cell problem that handles complex geometry in the interface in a natural fashion yet reduces to the well-known Ghost Fluid Method in simple cases. We test the accuracy of the estimates in a series of examples.

In biological cells, molecules are transported actively or by diffusion and react with each other when they are close.

The reactions occur with certain probability and there are few molecules of some chemical species. Therefore, a stochastic model is more accurate compared to a deterministic, macroscopic model for the concentrations based on partial differential equations.

At the mesoscopic level, the domain is partitioned into voxels or compartments. The molecules may react with other molecules in the same voxel and move between voxels by diffusion or active transport. At a finer, microscopic level, each individual molecule is tracked, it moves by Brownian motion and reacts with other molecules according to the Smoluchowski equation. The accuracy and efficiency of the simulations are improved by coupling the two levels and only using the micro model when it is necessary for the accuracy or when a meso description is unknown.

Algorithms for simulations with the mesoscopic, microscopic and meso-micro models will be described and applied to systems in molecular biology in three space dimensions.