Using argument notation to engineer biological simulations with increased confidence

The application of computational and mathematical modelling to explore the mechanics of biological systems is becoming prevalent. To significantly impact biological research, notably in developing novel therapeutics, it is critical that the model adequately represents the captured system. Confidence in adopting in silico approaches can be improved by applying a structured argumentation approach, alongside model development and results analysis. We propose an approach based on argumentation from safety-critical systems engineering, where a system is subjected to a stringent analysis of compliance against identified criteria. We show its use in examining the biological information upon which a model is based, identifying model strengths, highlighting areas requiring additional biological experimentation and providing documentation to support model publication. We demonstrate our use of structured argumentation in the development of a model of lymphoid tissue formation, specifically Peyer''s Patches. The argumentation structure is captured using Artoo (www.york.ac.uk/ycil/software/artoo), our Web-based tool for constructing fitness-for-purpose arguments, using a notation based on the safety-critical goal structuring notation. We show how argumentation helps in making the design and structured analysis of a model transparent, capturing the reasoning behind the inclusion or exclusion of each biological feature and recording assumptions, as well as pointing to evidence supporting model-derived conclusions.     

Role of cell polarity dynamics and motility in pattern formation due to contact-dependent signalling

A key challenge in biology is to understand how spatio-temporal patterns and structures arise during the development of an organism. An initial aggregate of spatially uniform cells develops and forms the differentiated structures of a fully developed organism. On the one hand, contact-dependent cell–cell signalling is responsible for generating a large number of complex, self-organized, spatial patterns in the distribution of the signalling molecules. On the other hand, the motility of cells coupled with their polarity can independently lead to collective motion patterns that depend on mechanical parameters influencing tissue deformation, such as cellular elasticity, cell–cell adhesion and active forces generated by actin and myosin dynamics. Although modelling efforts have, thus far, treated cell motility and cell–cell signalling separately, experiments in recent years suggest that these processes could be tightly coupled. Hence, in this paper, we study how the dynamics of cell polarity and migration influence the spatiotemporal patterning of signalling molecules. Such signalling interactions can occur only between cells that are in physical contact, either directly at the junctions of adjacent cells or through cellular protrusional contacts. We present a vertex model which accounts for contact-dependent signalling between adjacent cells and between non-adjacent neighbours through long protrusional contacts that occur along the orientation of cell polarization. We observe a rich variety of spatiotemporal patterns of signalling molecules that is influenced by polarity dynamics of the cells, relative strengths of adjacent and non-adjacent signalling interactions, range of polarized interaction, signalling activation threshold, relative time scales of signalling and polarity orientation, and cell motility. Though our results are developed in the context of Delta–Notch signalling, they are sufficiently general and can be extended to other contact dependent morpho-mechanical dynamics.     

Cost–benefit analysis of the mechanisms that enable migrating cells to sustain motility upon changes in matrix environments

Cells can move through extracellular environments with varying geometries and adhesive properties. Adaptation to these differences is achieved by switching between different modes of motility, including lamellipod-driven and blebbing motility. Further, cells can modulate their level of adhesion to the extracellular matrix (ECM) depending on both the level of force applied to the adhesions and cell intrinsic biochemical properties. We have constructed a computational model of cell motility to investigate how motile cells transition between extracellular environments with varying surface continuity, confinement and adhesion. Changes in migration strategy are an emergent property of cells as the ECM geometry and adhesion changes. The transition into confined environments with discontinuous ECM fibres is sufficient to induce shifts from lamellipod-based to blebbing motility, while changes in confinement alone within a continuous geometry are not. The geometry of the ECM facilitates plasticity, by inducing shifts where the cell has high marginal gain from a mode change, and conserving persistency where the cell can continue movement regardless of the motility mode. This regulation of cell motility is independent of global changes in cytoskeletal properties, but requires locally higher linkage between the actin network and the plasma membrane at the cell rear, and changes in internal cell pressure. In addition to matrix geometry, we consider how cells might transition between ECM of different adhesiveness. We find that this requires positive feedback between the forces cells apply on the adhesion points, and the strength of the cell–ECM adhesions on those sites. This positive feedback leads to the emergence of a small number of highly adhesive cores, similar to focal adhesions. While the range of ECM adhesion levels the cell can invade is expanded with this feedback mechanism; the velocities are lowered for conditions where the positive feedback is not vital. Thus, plasticity of cell motility sacrifices the benefits of specialization, for robustness.     

Model for self-polarization and motility of keratocyte fragments

Computational modelling of cell motility on substrates is a formidable challenge; regulatory pathways are intertwined and forces that influence cell motion are not fully quantified. Additional challenges arise from the need to describe a moving deformable cell boundary. Here, we present a simple mathematical model coupling cell shape dynamics, treated by the phase-field approach, to a vector field describing the mean orientation (polarization) of the actin filament network. The model successfully reproduces the primary phenomenology of cell motility: discontinuous onset of motion, diversity of cell shapes and shape oscillations. The results are in qualitative agreement with recent experiments on motility of keratocyte cells and cell fragments. The asymmetry of the shapes is captured to a large extent in this simple model, which may prove useful for the interpretation of experiments.     

Pattern selection and hysteresis in the Rietkerk model for banded vegetation in semi-arid environments

Banded vegetation is a characteristic feature of semi-arid environments. It occurs on gentle slopes, with alternating stripes of vegetation and bare ground running parallel to the contours. A number of mathematical models have been proposed to investigate the mechanisms underlying these patterns, and how they might be affected by changes in environmental conditions. One of the most widely used models is due to Rietkerk and co-workers, and is based on a water redistribution hypothesis, with the key feedback being that the rate of rainwater infiltration into the soil is an increasing function of plant biomass. Here, for the first time, we present a detailed study of the existence and stability of pattern solutions of the Rietkerk model on slopes, using the software package wavetrain (www.ma.hw.ac.uk/wavetrain). Specifically, we calculate the region of the rainfall–migration speed parameter plane in which patterns exist, and the sub-region in which these patterns are stable as solutions of the model partial differential equations. We then perform a detailed simulation-based study of the way in which patterns evolve when the rainfall parameter is slowly varied. This reveals complex behaviour, with sudden jumps in pattern wavelength, and hysteresis; we show that these jumps occur when the contours of constant pattern wavelength leave the parameter region giving stable patterns. Finally, we extend our results to the case in which a diffusion term for surface water is added to the model equations. The parameter regions for pattern existence and stability are relatively insensitive to small or moderate levels of surface water diffusion, but larger diffusion coefficients significantly change the subdivision into stable and unstable patterns.     

A parallel frontal solver for finite element applications

In finite element simulations, the overall computing time is dominated by the time needed to solve large sparse linear systems of equations. We report on the design and development of a parallel frontal code that can significantly reduce the wallclock time needed for the solution of these systems. The algorithm used is based on dividing the finite element domain into subdomains and applying the frontal method to each subdomain in parallel. The so‐called multiple front approach is shown to reduce the amount of work and memory required compared with the frontal method and, when run on a small number of processes, achieves good speedups. The code, HSL_MP42, has been developed for the Harwell Subroutine Library (http://www.numerical.rl.ac.uk/hsl). It is written in Fortran 90 and, by using MPI for message passing, achieves portability across a wide range of modern computer architectures. Copyright © 2001 John Wiley & Sons, Ltd.     

Application of the material force method to isotropic continuum damage

 The objective of this work is the exploitation of the notion of material forces in computational continuum damage mechanics. To this end we consider the framework of isotropic geometrically non–linear continuum damage and investigate the spatial and material settings that lead to either spatial or material forces, respectively. Thereby material forces essentially represent the tendency of material defects to move relative to the ambient material. In this work we combine an internal variable approach towards damage mechanics with the material force method. Thus the appearance of distributed material volume forces that are conjugated to the damage field necessitates the discretization of the damage variable as an independent field in addition to the deformation field. Consequently we propose a monolithic solution strategy for the corresponding coupled problem. The underlying kinematics, strong and weak forms of the coupled problem will be presented and implemented within a standard Galerkin finite element procedure. As a result in particular global discrete nodal quantities, the so–called material node point (surface) forces, are obtained and are studied for a number of computational examples. Received: 19 August 2002 / Accepted: 16 October 2002     

Importance of initial aortic properties on the evolving regional anisotropy, stiffness and wall thickness of human abdominal aortic aneurysms

Complementary advances in medical imaging, vascular biology and biomechanics promise to enable computational modelling of abdominal aortic aneurysms to play increasingly important roles in clinical decision processes. Using a finite-element-based growth and remodelling model of evolving aneurysm geometry and material properties, we show that regional variations in material anisotropy, stiffness and wall thickness should be expected to arise naturally and thus should be included in analyses of aneurysmal enlargement or wall stress. In addition, by initiating the model from best-fit material parameters estimated for non-aneurysmal aortas from different subjects, we show that the initial state of the aorta may influence strongly the subsequent rate of enlargement, wall thickness, mechanical behaviour and thus stress in the lesion. We submit, therefore, that clinically reliable modelling of the enlargement and overall rupture-potential of aneurysms may require both a better understanding of the mechanobiological processes that govern the evolution of these lesions and new methods of determining the patient-specific state of the pre-aneurysmal aorta (or correlation to currently unaffected portions thereof) through knowledge of demographics, comorbidities, lifestyle, genetics and future non-invasive or minimally invasive tests.     

On material forces and finite element discretizations

 The idea of using material forces also termed configurational forces in a computational setting is presented. The theory of material forces is briefly recast in the terms of a non-linear elastic solid. It is shown, how in a computational setting with finite elements (FE) the discrete configurational forces are calculated once the classical field equations are solved. This post-process calculation is performed in a way, which is consistent with the approximation of the classical field equations. Possible physical meanings of this configurational forces are discussed. A purely computational aspect of material forces is pointed out, where material forces act as an indicator to obtain softer discretizations. Received 12 December 2001 / Accepted 18 March 2002     

Modelling and analysis of bacterial tracks suggest an active reorientation mechanism in Rhodobacter sphaeroides

Most free-swimming bacteria move in approximately straight lines, interspersed with random reorientation phases. A key open question concerns varying mechanisms by which reorientation occurs. We combine mathematical modelling with analysis of a large tracking dataset to study the poorly understood reorientation mechanism in the monoflagellate species Rhodobacter sphaeroides. The flagellum on this species rotates counterclockwise to propel the bacterium, periodically ceasing rotation to enable reorientation. When rotation restarts the cell body usually points in a new direction. It has been assumed that the new direction is simply the result of Brownian rotation. We consider three variants of a self-propelled particle model of bacterial motility. The first considers rotational diffusion only, corresponding to a non-chemotactic mutant strain. Two further models incorporate stochastic reorientations, describing ‘run-and-tumble’ motility. We derive expressions for key summary statistics and simulate each model using a stochastic computational algorithm. We also discuss the effect of cell geometry on rotational diffusion. Working with a previously published tracking dataset, we compare predictions of the models with data on individual stopping events in R. sphaeroides. This provides strong evidence that this species undergoes some form of active reorientation rather than simple reorientation by Brownian rotation.     

Adaptive moving mesh methods for simulating one‐dimensional groundwater problems with sharp moving fronts

Accurate modelling of groundwater flow and transport with sharp moving fronts often involves high computational cost, when a fixed/uniform mesh is used. In this paper, we investigate the modelling of groundwater problems using a particular adaptive mesh method called the moving mesh partial differential equation approach. With this approach, the mesh is dynamically relocated through a partial differential equation to capture the evolving sharp fronts with a relatively small number of grid points. The mesh movement and physical system modelling are realized by solving the mesh movement and physical partial differential equations alternately. The method is applied to the modelling of a range of groundwater problems, including advection dominated chemical transport and reaction, non‐linear infiltration in soil, and the coupling of density dependent flow and transport. Numerical results demonstrate that sharp moving fronts can be accurately and efficiently captured by the moving mesh approach. Also addressed are important implementation strategies, e.g. the construction of the monitor function based on the interpolation error, control of mesh concentration, and two‐layer mesh movement. Copyright © 2002 John Wiley & Sons, Ltd.     

Simulation of surface pitting due to contact loading

A computational model for simulation of surface pitting of mechanical elements subjected to rolling and sliding contact conditions is presented. The two-dimensional computational model is restricted to modelling of high-precision mechanical components with fine surface finishing and good lubrication, where the cracks leading to pitting are initiated in the area of largest contact stresses at certain depth under the contacting surface. Hertz contact conditions with addition of friction forces are assumed and the position and magnitude of the maximum equivalent stress is determined by the finite element method. When the maximum equivalent stress exceeds the local material strength, it is assumed that the initial crack develops along the slip line in a single-crystal grain. The Virtual Crack Extension method in the framework of finite element analysis is then used for two-dimensional simulation of the fatigue crack propagation under contact loading from the initial crack up to the formation of the surface pit. The pit shapes and relationships between the stress intensity factor and crack length are determined for various combinations of contacting surface curvatures and loadings. The model is applied to simulation of surface pitting of two meshing gear teeth. Numerically predicted pit shapes in the face of gear teeth show a good agreement with the experimental observations. © 1998 John Wiley & Sons, Ltd.     

Dual finite element formulations for lumped reluctances coupling

A method for coupling magnetostatic and magnetodynamic finite element formulations with lumped reluctances is developed. Two dual h and b-conform formulations are extended to define the necessary magnetic relations between the magnetic fluxes and the magnetomotive forces. Adequate surface scalar potentials are defined and adequate boundary terms in the weak formulations are used. The magnetic relations can then be included in a reluctance network, in which the coupling of finite element regions and lumped regions aims at higher computational efficiency. The methods are developed in three dimensions, using a coupling of nodal and edge finite element approximations for the unknowns, and can easily be particularized in two dimensions.     

Application of fast Fourier transform cross-correlation for the alignment of large chromatographic and spectral datasets

Preprocessing of chromatographic and spectral data is an important aspect of analytical sciences. In particular, recent advances in proteomics have resulted in the generation of large data sets that require analysis. To assist accurate comparison of chemical signals, we propose two methods for the alignment of multiple spectral data sets. Based on methods previously described, each chromatograph or spectrum to be aligned is divided and aligned as individual segments to a reference. However, our methods make use of fast Fourier transform for the rapid computation of a cross-correlation function that enables alignments between samples to be optimized. The proposed methods are demonstrated in comparison with an existing method on a chromatographic and a mass spectral data set. It is shown that our methods provide an advantage of speed and a reduction of the number of input parameters required. The software implementations for the proposed alignment methods are available under the downloads section at http://ptcl.chem.ox.ac.uk/~jwong/specalign.     

Level set approach for the simulation of focused ion beam processing on the micro/nano scale

The level set method, introduced by Osher and Sethian (1988 J. Comput. Phys. 79 12-49), has recently become popular in the simulation of etching, deposition and photolithography processes in semiconductor manufacturing, as it is a highly robust and accurate computational technique for tracking moving interfaces. In this paper, the level set approach is applied to focused ion beam fabrication, allowing for the first time the simulation of targets with sub-regions that change their connectivity during processing. It is implemented in the code AMADEUS-level set (advanced modelling and design environment for sputter processes), which is capable of simulating surface topography changes in two dimensions taking re-deposition fluxes into account. We present two examples of comparisons between simulation and experiment that demonstrate the predictive capability of the code.     

Modelling platelet�Cblood flow interaction using the subcellular element Langevin method

In this paper, a new three-dimensional modelling approach is described for studying fluid–viscoelastic cell interaction, the subcellular element Langevin (SCEL) method, with cells modelled by subcellular elements (SCEs) and SCE cells coupled with fluid flow and substrate models by using the Langevin equation. It is demonstrated that: (i) the new method is computationally efficient, scaling as 𝒪(N) for N SCEs; (ii) cell geometry, stiffness and adhesivity can be modelled by directly relating parameters to experimentally measured values; (iii) modelling the fluid–platelet interface as a surface leads to a very good correlation with experimentally observed platelet flow interactions. Using this method, the three-dimensional motion of a viscoelastic platelet in a shear blood flow was simulated and compared with experiments on tracking platelets in a blood chamber. It is shown that the complex platelet-flipping dynamics under linear shear flows can be accurately recovered with the SCEL model when compared with the experiments. All experimental details and electronic supplementary material are archived at http://biomath.math.nd.edu/scelsupplementaryinformation/.     

Vacuum-ultraviolet reflectance difference spectroscopy for characterizing dielectrics-semiconductor interfaces

Reflectance difference spectroscopy (RDS) was applied to the characterization of SiO₂/Si, GeO₂/Ge, and high-k/III-V interface structures. We extended the spectral range of RDS to 8.4 eV in order to explore the optical transitions at the dielectrics-semiconductor interfaces as well as to have a high sensitivity to the interface anisotropy. Si surfaces with (110), (113), (331) and (120) orientations showed oxidation-induced RD changes in the vacuum-ultraviolet (VUV) range which were dependent on the surface orientation, oxidation method (dry or wet), and oxidation temperature. The Ge(110) surface also showed characteristic oxidation-induced changes in the VUV range, whereas Al₂O₃ deposition on GaAs(001) and InP(001) surfaces induced only the RD amplitude changes.     

Discontinuous modelling of masonry bridges

 Computational modelling frameworks for masonry bridges range from highly simplified methods to complex nonlinear finite element or discrete elements. In majority of cases the macro level nonlinear finite element models¹ and homogenisation techniques are adopted. Attention has also been given to assessment methodologies (discrete element method, rigid block spring method, lattice modelling, discontinuous deformation analysis, combined discrete/finite elements), which deal more directly with the discontinuous nature of structural masonry in a simplified micro modeling manner. These methods model an inherently discontinuous medium, but are also applied to problems where the transition from a continuum to discontinuum is important. Principal computational issue is the treatment of large number of distinct interacting domains, where the contact conditions are continuously updated and enforced as the solution progresses. Modelling of masonry arches requires a consideration of deformable multi-bodies and their contact nonlinearity, which is here realised in the context of the discontinuous deformation analysis, based on an assumed deformation field within distinct domains of arbitrary shapes with a rigorous imposition of contact constraints. Dedicated to the memory of Prof. Mike Crisfield, for his cheerfulness and cooperation as a colleague and friend over many years.     

Deciphering the flow structure of Czochralski melt using Partially Averaged Navier–Stokes (PANS) method

Czochralski melt flow is an outcome of complex interactions of centrifugal, buoyancy, coriolis and surface tension forces, which act at different length and time scales. As a consequence, the characteristic flow structures that develop in the melt are delineated in terms of recirculating flow cells typical of rotating Bénard–Marangoni convection. In the present study, Partially Averaged Navier–Stokes (PANS) method is used for the first time to study an idealized Czochralski crystal growth set-up. It is observed that with a reduction in the PANS filter width, more turbulent scales are resolved and the present PANS model is able to resolve almost all the characteristic flow structures in the Czochralski flow at a comparatively lower computational cost compared with more advanced turbulence modelling tools, such as Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES).