A new model for biological pattern formation

Various non-equilibrium growth models have been used to explore the development of morphology in biological systems. Here we review a class of biological growth models which exhibit fractal structures and discuss the relationship of these models to a variety of other phenomena.     

A finite element simulation scheme for biological muscular hydrostats

An explicit finite element scheme is developed for biological muscular hydrostats such as squid tentacles, octopus arms and elephant trunks. The scheme is implemented by embedding muscle fibers in finite elements. In any given element, the fiber orientation can be assigned arbitrarily and multiple muscle directions can be simulated. The mechanical stress in each muscle fiber is the sum of active and passive parts. The active stress is taken to be a function of activation state, muscle fiber shortening velocity and fiber strain; while the passive stress depends only on the strain. This scheme is tested by simulating extension of a squid tentacle during prey capture; our numerical predictions are in close correspondence with existing experimental results. It is shown that the present finite element scheme can successfully simulate more complex behaviors such as torsion of a squid tentacle and the bending behavior of octopus arms or elephant trunks.     

A theory of biological pattern formation

One of the elementary processes in morphogenesis is the formation of a spatial pattern of tissue structures, starting from almost homogeneous tissue. It will be shown that relatively simple molecular mechanisms based on auto- and cross catalysis can account for a primary pattern of morphogens to determine pattern formation of the tissue. The theory is based on short range activation, long range inhibition, and a distinction between activator and inhibitor concentrations on one hand, and the densities of their sources on the other. While source density is expected to change slowly, e.g. as an effect of cell differentiation, the concentration of activators and inhibitors can change rapidly to establish the primary pattern; this results from auto- and cross catalytic effects on the sources, spreading by diffusion or other mechanisms, and degradation.Employing an approximative equation, a criterium is derived for models, which lead to a striking pattern, starting from an even distribution of morphogens, and assuming a shallow source gradient. The polarity of the pattern depends on the direction of the source gradient, but can be rather independent of other features of source distribution. Models are proposed which explain size regulation (constant proportion of the parts of the pattern irrespective of total size). Depending on the choice of constants, aperiodic patterns, implying a one-to-one correlation between morphogen concentration and position in the tissue, or nearly periodic patterns can be obtained. The theory can be applied not only to multicellular tissues, but also to intracellular differentiation, e.g. of polar cells.The theory permits various molecular interpretations. One of the simplest models involves bimolecular activation and monomolecular inhibition. Source gradients may be substituted by, or added to, sink gradients, e.g. of degrading enzymes. Inhibitors can be substituted by substances required for, and depleted by activation.Sources may be either synthesizing systems or particulate structures releasing activators and inhibitors.Calculations by computer are presented to exemplify the main features of the theory proposed. The theory is applied to quantitative data on hydra — a suitable one-dimensional model for pattern formation — and is shown to account for activation and inhibition of secondary head formation.     

Anisotropic adaptive finite element method for modelling blood flow

In this study, we present an adaptive anisotropic finite element method (FEM) and demonstrate how computational efficiency can be increased when applying the method to the simulation of blood flow in the cardiovascular system. We use the SUPG formulation for the transient 3D incompressible Navier-Stokes equations which are discretised by linear finite elements for both the pressure and the velocity field. Given the pulsatile nature of the flow in blood vessels we have pursued adaptivity based on the average flow over a cardiac cycle. Error indicators are derived to define an anisotropic mesh metric field. Mesh modification algorithms are used to anisotropically adapt the mesh according to the desired size field. We demonstrate the efficiency of the method by first applying it to pulsatile flow in a straight cylindrical vessel and then to a porcine aorta with a stenosis bypassed by a graft. We demonstrate that the use of an anisotropic adaptive FEM can result in an order of magnitude reduction in computing time with no loss of accuracy compared to analyses obtained with uniform meshes.     

Statistical finite element method for real-time tissue mechanics analysis

The finite element (FE) method can accurately calculate tissue deformation. However, its low speed renders it ineffective for many biomedical applications involving real-time data processing. To accelerate FE analysis, we introduce a novel tissue mechanics simulation technique. This technique is suitable for real-time estimation of tissue deformation of specific organs, which is required in computer-aided diagnostic or therapeutic procedures. In this method, principal component analysis is used to describe each organ shape and its corresponding FE field for a pool of patients by a small number of weight factors. A mapping function is developed to relate the parameters of organ shape to their FE field counterpart. We show that irrespective of the complexity of the tissue's constitutive law or its loading conditions, the proposed technique is highly accurate and fast in estimating the FE field. Average deformation errors of less than 2% demonstrate the accuracy of the proposed technique.     

Anisotropic adaptive finite element method for modelling blood flow

In this study, we present an adaptive anisotropic finite element method (FEM) and demonstrate how computational efficiency can be increased when applying the method to the simulation of blood flow in the cardiovascular system. We use the SUPG formulation for the transient 3D incompressible Navier–Stokes equations which are discretised by linear finite elements for both the pressure and the velocity field.

Given the pulsatile nature of the flow in blood vessels we have pursued adaptivity based on the average flow over a cardiac cycle. Error indicators are derived to define an anisotropic mesh metric field. Mesh modification algorithms are used to anisotropically adapt the mesh according to the desired size field. We demonstrate the efficiency of the method by first applying it to pulsatile flow in a straight cylindrical vessel and then to a porcine aorta with a stenosis bypassed by a graft. We demonstrate that the use of an anisotropic adaptive FEM can result in an order of magnitude reduction in computing time with no loss of accuracy compared to analyses obtained with uniform meshes.     

Requirements for comparing the performance of finite element models of biological structures

The widespread availability of three-dimensional imaging and computational power has fostered a rapid increase in the number of biologists using finite element analysis (FEA) to investigate the mechanical function of living and extinct organisms. The inevitable rise of studies that compare finite element models brings to the fore two critical questions about how such comparative analyses can and should be conducted: (1) what metrics are appropriate for assessing the performance of biological structures using finite element modeling? and, (2) how can performance be compared such that the effects of size and shape are disentangled? With respect to performance, we argue that energy efficiency is a reasonable optimality criterion for biological structures and we show that the total strain energy (a measure of work expended deforming a structure) is a robust metric for comparing the mechanical efficiency of structures modeled with finite elements. Results of finite element analyses can be interpreted with confidence when model input parameters (muscle forces, detailed material properties) and/or output parameters (reaction forces, strains) are well-documented by studies of living animals. However, many researchers wish to compare species for which these input and validation data are difficult or impossible to acquire. In these cases, researchers can still compare the performance of structures that differ in shape if variation in size is controlled. We offer a theoretical framework and empirical data demonstrating that scaling finite element models to equal force: surface area ratios removes the effects of model size and provides a comparison of stress-strength performance based solely on shape. Further, models scaled to have equal applied force:volume ratios provide the basis for strain energy comparison. Thus, although finite element analyses of biological structures should be validated experimentally whenever possible, this study demonstrates that the relative performance of un-validated models can be compared so long as they are scaled properly.     

A one-dimensional finite element method for simulation-based medical planning for cardiovascular disease

We have previously described a new approach to planning treatments for cardiovascular disease, Simulation-Based Medical Planning, whereby a physician utilizes computational tools to construct and evaluate a combined anatomic/physiologic model to predict the outcome of alternative treatment plans for an individual patient. Current systems for Simulation-Based Medical Planning utilize finite element methods to solve the time-dependent, three-dimensional equations governing blood flow and provide detailed data on blood flow distribution, pressure gradients and locations of flow recirculation, low wall shear stress and high particle residence. However, these methods are computationally expensive and often require hours of time on parallel computers. This level of computation is necessary for obtaining detailed information about blood flow, but likely is unnecessary for obtaining information about mean flow rates and pressure losses. We describe, herein, a space-time finite element method for solving the one-dimensional equations of blood flow. This method is applied to compute flow rate and pressure in a single segment model, a bifurcation, an idealized model of the abdominal aorta, in three alternate treatment plans for a case of aorto-iliac occlusive disease and in a vascular bypass graft. All of these solutions were obtained in less than 5 min of computation time on a personal computer.     

A penetration-based finite element method for hyperelastic 3D biphasic tissues in contact. Part II: finite element simulations

The penetration method allows for the efficient finite element simulation of contact between soft hydrated biphasic tissues in diarthrodial joints. Efficiency of the method is achieved by separating the intrinsically nonlinear contact problem into a pair of linked biphasic finite element analyses, in which an approximate, spatially and temporally varying contact traction is applied to each of the contacting tissues. In Part I of this study, we extended the penetration method to contact involving nonlinear biphasic tissue layers, and demonstrated how to derive the approximate contact traction boundary conditions. The traction derivation involves time and space dependent natural boundary conditions, and requires special numerical treatment. This paper (Part II) describes how we obtain an efficient nonlinear finite element procedure to solve for the biphasic response of the individual contacting layers. In particular, alternate linearization of the nonlinear weak form, as well as both velocity-pressure, v-p, and displacement-pressure, u-p, mixed formulations are considered. We conclude that the u-p approach, with linearization of both the material law and the deformation gradients, performs best for the problem at hand. The nonlinear biphasic contact solution will be demonstrated for the motion of the glenohumeral joint of the human shoulder joint.     

A new method for finite element simulation of orthodontic appliance-teeth-periodontium-alveolus system

This paper describes a new simulation method to analyze the initial behavior of the total system comprising orthodontic appliance, teeth, and their supporting structures. It is based on a finite element method which additionally takes account of a rotational degree of freedom. Beam and rod elements are used for finite element idealization of orthodontic appliance. Through spring elements it is connected with the teeth supported by the alveolar structures. The technique of 'initial strain' is introduced so as to analyze the effects of a gable bend and activation on the force system which is delivered by the orthodontic appliance. As compared with the photoelastic technique hitherto used, this method serves to investigate systematically and quantitatively the initial aspect of orthodontic tooth movement.     

Tissue-fluid interface analysis using biphasic finite element method

Numerical studies on fluid-structure interaction have primarily relied on decoupling the solid and fluid sub-domains with the interactions treated as external boundary conditions on the individual sub-domains. The finite element applications for the fluid-structure interactions can be divided into iterative algorithms and sequential algorithms. In this paper, a new computational methodology for the analysis of tissue-fluid interaction problems is presented. The whole computational domain is treated as a single biphasic continuum, and the same space and time discretisation is carried out for the sub-domains using a penalty-based finite element model. This procedure does not require the explicit modelling of additional boundary conditions or interface elements. The developed biphasic interface finite element model is used in analysing blood flow through normal and stenotic arteries. The increase in fluid flow velocity when passing through a stenosed artery and the drop in pressure at the region are captured using this method.     

Effects of grid dimensions on finite element models of an articular surface

Recent activity in finite element analysis of articular joints has emphasized refinements in geometry and material properties. In the implementation of such models, it is necessary to ensure that grid dimensions are optimal for suitable solutions of displacements, strains and stresses. A method of grid optimization was developed to ensure that for typical material properties, finite element models of an articular surface agree with known analytical solutions. The layered axisymmetric model presented by Askew and Mow (J. biomech. Engng 100, 105-115, 1978) was used as a reference. From this reference, an STZ of 0.2 mm, middle and deep zones of 0.8 mm and tidemark region of 0.2 mm were chosen. Cancellous bone was an infinite elastic half space under these layers. Loading was a parabolic distribution over a 10 mm radius having a peak of 1 MPa. Agreement was obtained between analytical solutions and finite element solutions when the finite element model had a radial boundary of 30 mm radius and a bone depth of 32 mm. These results suggested that in models of real joints, care must be taken to ensure the boundaries are reasonably represented and that sufficient bone is modelled for adequate solutions.     

Using mathematical models to help understand biological pattern formation

One of the characteristics of biological systems is their ability to produce and sustain spatial and spatio-temporal pattern. Elucidating the underlying mechanisms responsible for this phenomenon has been the goal of much experimental and theoretical research. This paper illustrates this area of research by presenting some of the mathematical models that have been proposed to account for pattern formation in biology and considering their implications.     

On the stationary state analysis of reaction-diffusion mechanisms for biological pattern formation

We present a biologically plausible two-variable reaction-diffusion model for the developing vertebrate limb, for which we postulate the existence of a stationary solution. A consequence of this assumption is that the stationary state depends on only a single concentration-variable. Under these circumstances, features of potential biological significance, such as the dependence of the steady-state concentration profile of this variable on parameters such as tissue size and shape, can be studied without detailed information about the rate functions. As the existence and stability of stationary solutions, which must be assumed for any biochemical system governing morphogenesis, cannot be investigated without such information, an analysis is made of the minimal requirements for stable, stationary non-uniform solutions in a general class of reaction-diffusion systems. We discuss the strategy of studying stationary-state properties of systems that are incompletely specified. Where abrupt transitions between successive compartment-sizes occur, as in the developing limb, we argue that it is reasonable to model pattern reorganization as a sequence of independent stationary states.     

Dynamic finite element implementation of nonlinear,anisotropic hyperelastic biological membranes

We present a novel method for the implementation of hyperelastic finite strain, non-linear strain-energy functions for biological membranes in an explicit finite element environment. The technique is implemented in LS-DYNA but may also be implemented in any suitable non-linear explicit code. The constitutive equations are implemented on the foundation of a co-rotational uniformly reduced Hughes-Liu shell. This shell is based on an updated-Lagrangian formulation suitable for relating Cauchy stress to the rate-of-deformation, i.e. hypo-elasticity. To accommodate finite deformation hyper-elastic formulations, a co-rotational deformation gradient is assembled over time, resulting in a formulation suitable for pseudo-hyperelastic constitutive equations that are standard assumptions in biomechanics. Our method was validated by comparison with (1) an analytic solution to a spherically-symmetric dynamic membrane inflation problem, incorporating a Mooney-Rivlin hyperelastic equation and (2) with previously published finite element solutions to a non-linear transversely isotropic inflation problem. Finally, we implemented a transversely isotropic strain-energy function for mitral valve tissue. The method is simple and accurate and is believed to be generally useful for anyone who wishes to model biologic membranes with an experimentally driven strain-energy function.     

Trabecular surface remodeling simulation for cancellous bone using microstructural voxel finite element models

A computational simulation method for three-dimensional trabecular surface remodeling was proposed, using voxel finite element models of cancellous bone, and was applied to the experimental data. In the simulation, the trabecular microstructure was modeled based on digital images, and its morphological changes due to surface movement at the trabecular level were directly expressed by removing/adding the voxel elements from/to the trabecular surface. A remodeling simulation at the single trabecular level under uniaxial compressive loading demonstrated smooth morphological changes even though the trabeculae were modeled with discrete voxel elements. Moreover, the trabecular axis rotated toward the loading direction with increasing stiffness, simulating functional adaptation to the applied load. In the remodeling simulation at the trabecular structural level, a cancellous bone cube was modeled using a digital image obtained by microcomputed tomography (microCT), and was uniaxially compressed. As a result, the apparent stiffness against the applied load increased by remodeling, in which the trabeculae reoriented to the loading direction. In addition, changes in the structural indices of the trabecular architecture coincided qualitatively with previously published experimental observations. Through these studies, it was demonstrated that the newly proposed voxel simulation technique enables us to simulate the trabecular surface remodeling and to compare the results obtained using this technique with the in vivo experimental data in the investigation of the adaptive bone remodeling phenomenon.     

A finite element method for mechanical response of hair cell ciliary bundles

This paper describes the development of a methodology for performing a mechanical analysis of hair cell ciliary bundles. The cilia were modeled as shear deformable beams, and interconnections were modeled as two-force members. These models were incorporated into software, which performs a finite element analysis of a user-defined bundle. The algorithm incorporates aspects of the bundle such as geometric realignment and buckling of compressed side links. A sample bundle is introduced and results of modeling it are presented.     

A fully implicit finite element method for bidomain models of cardiac electrophysiology

This work introduces a novel, unconditionally stable and fully coupled finite element method for the bidomain system of equations of cardiac electrophysiology. The transmembrane potential Φ(i)-Φ(e) and the extracellular potential Φ(e) are treated as independent variables. To this end, the respective reaction-diffusion equations are recast into weak forms via a conventional isoparametric Galerkin approach. The resultant nonlinear set of residual equations is consistently linearised. The method results in a symmetric set of equations, which reduces the computational time significantly compared to the conventional solution algorithms. The proposed method is inherently modular and can be combined with phenomenological or ionic models across the cell membrane. The efficiency of the method and the comparison of its computational cost with respect to the simplified monodomain models are demonstrated through representative numerical examples.