Introductory Slides YouTube Video Tutorial

Introduction to Multi-Scale Modeling in Continuity

Continuity is a problem solving environment for multi-scale modeling in biomechanics, biotransport and electrophysiology. It is freely available for academic use under license from UCSD and downloadable at the Continuity download page. Several new updates are released per year. Although Continuity is generic and application-neutral, we have developed it to support our research on cardiac physiology and pathophysiology, and therefore it is particularly well adapted to multi-scale problems in cardiac electrophysiology and biomechanics including coupled, multi-physics, multi-scale models of the heart. In this introductory guided tour, we focus on example applications in cardiac modeling.

Because multi-scale biological models are most often valuable when used in close association with experimental investigation, Continuity incorporates extensive facilities for analyzing and fitting measured data, including biomedical images, and the ability to author and dynamically compile component models, such as constitutive equations for tissue properties, systems models for cellular dynamics e.g. ionic current models in electrophysiology, material coordinate definitions, and lumped parameter models of circulatory hemodynamics. A key objective and design priority of our continuing Continuity development efforts is that every aspect of a problem or simulation that you would consider to constitute the mathematical model or its biological parameters (as distinct from the numerical implementation) should be data that can be easily shared and re-used. For that reason too, there is also a public Continuity model database that you can access directly from within the Continuity client. We are close to realizing that goal, though certain parts of the model, most notably the underlying partial differential equations (PDEs) themselves, are still currently expressed in Continuity source code.

At its core, Continuity is a PDE solver that implements finite element methods for nonlinear mechanics, reaction-diffusion and monodomain problems. A typical multi-scale problem in Continuity spans three or more scales, typically comprising a system of ODEs at each point in the physical domain representing a network of biochemical or biophysical interactions, a constitutive law representing physical properties of the system that spatially couple these pointwise networks, a system of PDEs representing the conservation law governing the physics over the anatomy of the system, and a model of boundary conditions representing the interactions of the system with neighboring structures. For example, in a typical cardiac biomechanics problem, ODEs represent dynamic Markov models of myofilament activation and interactions in the cardiac myocyte, the constitutive equations describe the three-dimensional anisotropic mechanical properties of the tissue with respect to local structural axes and how the myofilament interactions modify them dynamically, the governing PDEs for force and momentum balance are solved for an anatomic model of the heart subject to boundary conditions developed by a lumped-parameter model of the whole circulation. Similarly, a cardiac action potential propagation model comprises equations representing the dynamics of individual ion channels and whole cell ionic currents and potentials as systems of ODEs, that are spatially coupled via the monodomain PDEs by anisotropic conductivities and subject to boundary conditions representing the passive current spread through the torso.

Installing Continuity

Continuity can be downloaded for Windows, (Intel) Mac OS X and Linux platforms from the Continuity download page. Versions 7113 and later are preconfigured to allow the authoring and compilation of user-defined models such as material coordinate transformations and constitutive laws. Older versions require the Compiler Toolkit Plug-In for Windows or the Set-Up Build Environment Script for Mac OS X found here. The entire history of releases since Fall 2005 is there in case you need to go back to an earlier working version in the event of a problem in a subsequent release. There is also a subversion (svn) source code repository for developers. Installation uses industry-standard mechanisms specific to each platform, such as dmg files for the MacOS and an msi point and click installer for Windows. See the help pages for detailed installation instructions. When you run Continuity, you will see the option to register yourself as a user from the splash screen. You don’t need to register to install or use Continuity, but access to the model database does require a registration key as it makes use of federally supported Internet resources at UCSD and the San Diego Supercomputer Center. We encourage you to register and provide us with information about your research interests so that we can improve the tools for your needs and provide useful statistical information to our federal sponsors.

The Continuity Architecture

Continuity consists of a client and a server. The client serves as the graphical user interface, and the server carries out the numerical calculations for simulation and problem solving. The installation process installs both on your computer. When you launch Continuity the graphical user interface (GUI) client starts and gives you the default option of running the client and server together on your local machine so that they behave as a single application. However, you can also use the Continuity client to connect to a remote server running on another host including a Linux cluster. There is no need for the client and server to be running on the same OS platform. The server has dynamically compilable and loadable numerical analysis functions for mesh handling, data fitting and image handling, biomechanics, electrophysiology, and transport processes. Correspondingly, the user can load client-side components corresponding to these separate modules that appear as new menus in the GUI client. Indeed, you will see the option to load certain modules directly in the splash screen. By default only the Mesh module is selected. You don’t have to decide at this stage which module you want to use. You can always load modules at any time. For a more complete summary of Continuity’s design and architecture, see the Continuity overview.

Using the Mesh Module to Make a Simple Ventricular Model

The Mesh module is the most commonly used module in Continuity, which is why it is active by default when you start Continuity. (Note that you can use the Preferences menu to change the default modules that are selected when Continuity is launched). One powerful and uncommon feature of Continuity is the ability to create and use meshes using curvilinear coordinate systems, such as cylindrical polar or spherical polar coordinates, as an alternative to rectangular Cartesian coordinates. You will soon be able to define your own curvilinear world coordinate system in Continuity, but for the meantime, there are four built-in world coordinate systems including Cartesian coordinates.

Prolate Spheroidal Approximations of Ventricular Geometry

Investigators working on ventricular mechanics in the 1960’s and 70’s, including Streeter, Hanna and Mirsky, approximated the left ventricle as nested ellipses of revolution. Here we show data from Streeter and Hanna (1973), who fitted truncated ellipses of revolution to measurements of left ventricular epicardial and endocardial major and minor axis dimensions in canine hearts, isolated and fixed under various loading conditions.

Just as using cylindrical coordinates simplifies the description of tubular objects, the natural coordinate system for ellipsoidal shapes is prolate spheroidal (elliptic-hyperbolic-polar) coordinates (λ,μ,θ). In prolate spheroidal coordinates, the longitudinal (μ) coordinate forms families of ellipses, each with a common focal point at x=d. The radial (λ) coordinates form hyperbolae. Streeter and Hanna showed that, the focal distance is remarkably constant both between the epicardial and endocardial surfaces and between systole and diastole in the dog left ventricle. Based on their data we will create a prolate spheroidal mesh in Continuity using a focal length of 3.7 cm. From within the Continuity Mesh menu, choose Edit→Coordinates…, select Prolate Spheroidal in the pop-up menu, enter 3.7 for the focus position and click OK to submit.

• Continuity Command Summary From within the Mesh menu, choose Edit→Coordinates… Select Prolate Spheroidal in the pop-up menu Enter 3.7 for the focus position Click OK to submit.

   

From the major and minor radii reported by Streeter and Hanna, the focus, and the equations for prolate spheroidal coordinates we can easily compute the λ coordinates of the epicardial and endocardial surfaces, as shown in the table. Streeter and Hanna also defined a truncation factor as the fraction of the equator-apex distance between the base and the equator. Their estimate of nearly 0.5 translates to a μ coordinate at the base of 120°.

Finite Element Interpolation

We’ll use these values to specify nodes of our mesh, but first we need to choose a suitable finite element basis function to interpolate the geometric variables. Unlike finite difference methods, which discretize a continuum into a finite number of points and approximate the governing differential equations using differences, finite element methods approximate the continuum using a finite number of functions, usually piecewise polynomials spanning subdomains (elements) with parameters defined at the vertices (nodes). The governing equations are represented using a weak, integral form that gives rise to sums instead of differences. This facilitates the finite element assembly of the global system of equations by addition of the element matrices representing the discretization of the integral equations over each subdomain.

Continuity supports a variety of basis functions for isoparametric finite element interpolation. The simplest type of elements are piecewise linear Lagrange elements. Note that within each element, there is a local finite element (ξ) coordinate that always varies from 0 at the first node to 1 at the second. The parameters of the linear interpolation are the values (u₁ and u₂) of the function at the nodes, and thus the basis functions (φ₁ and φ₂)are the corresponding linear functions that when multiplied by u₁ and u₂, respectively, and added together result in the linear interpolation between u₁ and u₂. Note also that the term isoparametric here refers to the fact that rather than interpolate a variable y as a function of geometric variables x, we interpolate both as functions of ξ, so that y(x) does not itself need to be a function. The major advantage of finite element basis functions such as linear Lagrange interpolation schemes is that, since the parameters of the interpolation are the value at global nodes shared by adjacent elements, (C₀) continuity of interpolated variables between neighboring elements is automatically ensured.

Simply by taking tensor products, we can use these simple 1-D linear Lagrange basis functions to generate 2-D quadralateral bilinear Lagrange elements as functions of ξ₁ and ξ₂ and 3-D hexahedral trilinear Lagrange elements as functions of ξ₁, ξ₂, and ξ₃.

We are now ready to make our first mesh in Continuity, using the prolate spheroidal coordinates we calculated from Streeter and Hanna’s data. We will build a simple truncated, thick-walled prolate spheroidal mesh with a single linear 3-D element. First, we need to select a trilinear Lagrange 3-D basis function. Using the Edit→Basis… command from the Mesh menu, select Lagrange Basis Function→3D→Linear-Linear-Linear with 3 integration/collocation points for “Xi 1” (i.e. ξ₁), “Xi 2” (i.e. ξ₂), and “Xi 3” (i.e. ξ₃). Click Add and OK to submit.

• Continuity Command Summary Using the Edit→Basis… command from the Mesh menu, select Lagrange Basis Function→3D→Linear-Linear-Linear Use 3 integration/collocation points for “Xi 1” (i.e. ξ1), “Xi 2” (i.e. ξ2), and “Xi 3” (i.e. ξ3). Click Add then OK to submit.

   

Nodes

Selecting Edit→Nodes… command from the Mesh menu brings up a large form where we can define the nodal values of geometric coordinates and other pre-defined variables such as fiber angles and fields. In our model, Coordinates 1, 2 and 3 are λ, μ and θ, respectively. Even though a trilinear element has eight element nodes we only need four nodes for this model, because it is axisymmetric, wrapping around on itself. We’ll see how to handle this when we define elements, shortly. Since we only selected one basis function, Linear-Linear-Linear Lagrange 3*3*3 probably already appears at the top of each coordinate column in the Nodes form. If not, select it from the menu for all three coordinates. Now, click the Insert Node button on the left three times to so that there are a total of 4 nodes. Then simply type in the nodal prolate spheroidal coordinates for each node that we computed for the canine left ventricle. Note that the default units for angular coordinates such as μ and θ is degrees, but that can be changed to radians using the radio button at the bottom of the form. We will also assign a label to each node, which will be useful later when we define boundary conditions. To do this enter a text string in the text entry field at the top of the form next to the Edit button to the right of the “Node Number:” display. Use the labels APEX_ENDO, BASE_ENDO, APEX_EPI, BASE_EPI for the nodes 1-4, respectively, as shown in the table. The underscore in the label serves the special purpose of enabling a compound label, which we will take advantage of later when we refine the mesh. The grayed out fields next to the unchecked boxes for each nodal parameter are not needed for linear elements. They represent derivative parameters that we will use later for higher-order Hermite interpolation schemes. Special note: When you are entering nodal labels, you must press return for the entry to be saved. The Import/Export/Graph button on the left panel provides a convenient way to import or export the nodes as tab-delimited tables that can be viewed or edited in the Continuity Table Manager or a spreadsheet program, such as Microsoft Excel. Click OK to submit the nodal coordinates and labels.

• Continuity Command Summary Using the Edit→Nodes… command from the Mesh menu Make sure Linear-Linear-Linear Lagrange 3*3*3 appears in top of each coordinate Insert 3 nodes, so there are a total of 4 Define each node with values and labels Click OK to submit the nodal coordinates and labels

   

Elements

Next, we’ll use the Edit→Elements… command from the Mesh menu to create one element. By entering the eight node numbers for each trilinear element, we both define the elements and the directions of their local ξ coordinates. The graphic in the form helps remind us of the numbering convention in which nodes are entered starting at (ξ₁,ξ₂,ξ₃)=(0,0,0) and ending with the node at (ξ₁,ξ₂,ξ₃)=(1,1,1). By convention for meshes defined in curvilinear coordinates, we make ξ₁ coincide (to the extent possible) with the polar (θ) axis, ξ₂ with the “longitudinal” (z, φ or μ) axis, and ξ₃ with the “radial” (r or λ) axis, while also trying to keep the ξ coordinates right-handed by convention. Note that this has the odd effect of requiring that ξ₁ coincide with the negative θ axis in prolate spheroidal coordinates because θ is the third prolate spheroidal coordinate by convention, but it is the second coordinate in cylindrical or spherical polar coordinates. However, in this example, ξ₁ wraps all the way around, so the order doesn’t matter. Hence following these conventions, the eight element nodes of our single trilinear element are defined by the four global nodes (1-4) entering the node sequence 1,1,2,2,3,3,4,4 in the element form for element 1, where we have incremented first circumferentially (ξ₁), then longitudinally (ξ₂) from apex to base, and finally radially (ξ₃) from endocardium to epicardium. Note that in this very simple case, our world curvilinear (θ_j) coordinates and our local finite element (ξ_i) coordinates coincide and both are right-handed systems, though their coordinate order is different. In general, though, we will not be dealing with such simple shapes so we can not expect the ξ coordinates to follow along world coordinate axes, though we do generally try to maintain righted-handedness. You should also enter a label for Element 1 (e.g. “LV”) into the Element Label field at the top of the form (remember to press return). Click OK to submit the elements.

• Continuity Command Summary Using the Edit→Elements… command from the Mesh menu Enter the node sequence 1,1,2,2,3,3,4,4 in the element form for element 1 Enter a label for Element 1 (e.g. “LV”) into the Element Label field at the top of the form Click OK to submit the nodal coordinates and labels File→Save→Model to save your model

   

Calculating Mesh

Now we are ready to use our new mesh. The button on the toolbar is a shortcut to the Mesh→Calculate Mesh… menu choice. Click OK in the dialog box. Behind the scenes, Continuity also carried out a File→Send to Server command, which updates the Continuity server with any client changes and can be performed using the toolbar button. Nothing visible happens yet, but now you can render the mesh as a wire frame or surfaces using the Mesh→Render→Elements… command or the toolbar button. There is also a Mesh→Render→Nodes… command. There is a wide variety of graphical rendering controls in the View menu. For example, the View→Show→OpenMesh command lists the rendered objects and their attributes. The entry field next to “select” in the Information tab provides a way to highlight specific elements or nodes. Select the nodes object from the list of graphical objects on the left to try this.

• Continuity Command Summary Calculate the mesh, using either or Mesh→Calculate Mesh… Render elements using either or Mesh→Render→Elements…

   

Mesh Volumes

The Mesh→List→Elements… command will give you information on element volumes, faces, and lines. Check that the volume of the wall is being integrated accurately. If necessary, you can adjust the number of Gaussian quadrature integration points in the Edit→Basis… form to increase the accuracy of numerical integration, at the expense of more calculations per element. Although the left ventricular cavity itself is not part of the mesh, there are also ways to integrate the cavity volume bounded by defined surfaces in Continuity. Choose the Mesh→List→Volume… command. This command calculates volumes bounded by a surface surrounding a single origin. To get an accurate cavity volume, we need to place the origin at the plane of the endocardial base. We can use the information in the geometric data table to compute this position. Y and Z are zero but X is negative and equal to hxa, which is 0.45×4.2=1.89 for the endocardium at end-diastole. For (X,Y,Z) under “Define the Origin” enter (-1.89,0.0,0.0). Use the pull-down menu next to the “Element List” field and you should see your previously defined element label LV. Select it. Choose LV Endocardium from the menu underneath “Wall”, then click the Insert Surfaces button. Then select Epicardium from the “Wall” menu and Insert Surfaces again. Click OK and the cavity and wall volumes will be displayed. Because the calculation is completely different, the wall volume will probably not be exactly the same as the mesh volume displayed previously, but it should be fairly close.

• Continuity Command Summary Mesh→List→Elements… Mesh→List→Volume… “Define the Origin” enter (-1.89,0.0,0.0) Select LV from “Element List” field Select LV Endocardium from “Wall” field Click the Insert Surfaces Select Epicardium from the “Wall” Click the Insert Surfaces again Click OK

   

Ventricular Myofiber Architecture

The ventricular muscle fiber architecture has a characteristic pattern in most mammals that is surprisingly conserved. Roughly, the myofibers remain in planes almost parallel to the walls with helical pitch orientations varying from 60° to 90° counterclockwise from the circumferential axis (right-handed helices) in the subendocardium to -45° to -90° (left handed) in the epicardium. This familiar convention for fiber orientation is due to the quantitative measurements of Streeter and co-workers in the 1960’s, and for this fiber orientations reported using these angles are sometimes called Streeter angles. Most models assume that fibers lie parallel to the wall, though some investigators have measured small but significant out-of-plane transverse or imbrication angles, especially near the base and apex, which we will assume to be zero. It is important to remember that cardiac myocytes, unlike skeletal myocytes, are not really fibers at all; they are rod-shaped cells with end-to-end and side-to-side junctional connections and branches. Sections cut across the fibers reveal that the myofibers are organized by the perimysial collagen matrix into branching laminar myocardial sheets approximately four cells thick. Their orientations tend to have a greater local dispersion and commonly show bimodal distributions in some regions, but based on the predominant sheet orientations, it is now common in biomechanical and electrical models of the ventricles to use a local, orthonormal fiber-sheet-normal coordinate frame following the work of LeGrice and colleagues.

We will use published measurements of fiber and sheet angles to build the material coordinates in our model. These measurements show that muscle fiber and sheet orientations vary transmurally in the dog. A comprehensive analysis of fiber and sheet architecture in the canine ventricles was performed by LeGrice et al (1997). But even with our single element prolate spheroidal model, we can represent much of the variation in myofiber architecture. To do so we will need basis functions that are higher order than linear.

Quadratic Lagrange and Cubic Hermite Basis Functions

By placing a third node at the middle of each element (at ξ=0.5), we now have three parameters that we can use to interpolate second degree polynomials using the quadratic Lagrange basis functions. Note that each basis function is itself a quadratic having the property that it is one at the node corresponding to its own nodal parameter and zero at all of the other nodes. In the same way we can use four equally spaced nodes for cubic Lagrange interpolation, but a more useful strategy that requires only end (or vertex) nodes (at ξ=0 and 1) is to use the cubic Hermite basis functions. In 1-D, these four basis functions correspond to two parameters at each of the two end nodes, one being the value of the interpolated variable and the second being its first derivative with respect to ξ.

A big advantage of cubic Hermite interpolation is not only that the elements can now support a higher-order variation of mesh variables, but that now, by sharing global nodal parameters, we can now achieve first derivative (C₁) continuity of mesh variables between elements. A practical limitation arises here though because the ξ coordinates themselves are not generally continuous. For example, if neighboring elements are different in length (which is a common and useful property of typical unstructured finite element meshes), then shared nodal derivatives with respect to ξ will not result in smooth slopes in Euclidean space. We correct this problem in Continuity by storing derivatives with respect to global arc lengths s at global nodes, and converting them into local ξ derivatives in each element using local element scaling factors, which are computed from the element arc lengths when you issue the Mesh→Calculate Mesh… command in Continuity.

In the same manner as we did for linear Lagrange basis functions, we can construct 2-D and 3-D quadratic Lagrange and cubic Hermite elements by taking tensor products. We can even combine different 1-D basis functions in each direction, to create, e.g. a 6-noded quadratic-linear Lagrange 2-D element. Forming a 3-D tricubic Hermite element results in an 8-noded hexahedral element with 64 nodal parameters, comprising 8 derivative parameters per node. For the next step, we will need two new basis functions: a linear-bicubic tensor product basis function with 4 degrees of freedom per node, that we will use to interpolate transverse angles, and a bilinear-cubic basis function with two degrees of freedom per node for the sheet angles. Using the Edit→Basis… command from the Mesh menu, select Hermite Basis Function→3D→Linear-Linear-Cubic with 3 integration/collocation points for each local Xi coordinate direction as before. Click Add and OK to submit.

• Continuity Command Summary Using the Edit→Basis… command from the Mesh menu, select Hermite Basis Function→3D→Linear-Linear-Cubic Use 3 integration/collocation points for each local Xi coordinate. Click Add then OK to submit.

   

Fiber/Sheet Angles

For the purposes of demonstration here, we will approximate the measured transmural distributions of fiber and sheet angles in our 3-D prolate spheroidal model by using the trilinear Lagrange basis functions that you created. Later we will learn about the data Fitting module of Continuity which enables us to fit nodal parameters of meshes directly to measurements of geometry or field variables. We will use a trilinear basis function to interpolate fiber angles and estimate nodal values from the plotted data. Since they are more variable than fibers we will only specify values for the sheet angles, which we will also approximate using trilinear Lagrange interpolation. Estimated parameters are tabulated here. To enter these values, let’s go back to the Edit→Nodes… form, which you select from the Mesh menu. Click on the Fiber Angles tab. Then select Linear-Linear-Linear 3*3*3 from the Select Basis Number drop-down menu below “Fiber Angle”, select Linear-Linear-Linear 3*3*3 under “Transverse Angle”, and Linear-Linear-Linear 3*3*3 under “Sheet Angle”. Now you will see that extra nodal parameter entry fields are checked. Use the table to enter the nodal value parameters (ignore s(3)).

• Continuity Command Summary Mesh→Edit→Nodes… Click on the Fiber Angles tab Select Linear-Linear-Linear 3*3*3 under Fiber Angle Select Linear-Linear-Linear 3*3*3 under Transverse and Sheet Angle Enter values Click OK Save model

   

Material Coordinates

Continuity uses four coordinate systems: Global Rectangular Cartesian coordinates (Y₁,Y₂,Y₃); curvilinear World coordinates (Θ₁,Θ₂,Θ₃₁); local finite element coordinates (ξ₁,ξ₂,ξ₃); and local material coordinates (x₁,x₂,x₃). The Mesh→Edit Material Coordinates command allows us to define an orthogonal transformation that rotates global Cartesian coordinates at any point in the mesh to the local orthonormal material coordinates. The standard predefined material coordinate transformation in Continuity uses the local ξ coordinates as intermediates to create a locally orthogonal body coordinate system with one axis aligned with ξ₁, one in the ξ₁-ξ₂ plane and the third perpendicular to the other two. These axes are then rotated through the fiber, transverse and sheet angles to generate material coordinates coinciding with the fiber, sheet and sheet-normal directions. The fiber, transverse and sheet angles are the local varying parameters of this model. They carry no intrinsic meaning to Continuity and are only labeled as such in the Nodes editor for historical reasons when this material coordinate transformation was built-in and hard-coded in Continuity. Within the context of Continuity, we call this definition of material axis transformation the Standard Material Coordinates model. The equations and parameters of this material model have already been defined and are stored in the Continuity database.

Go to Mesh->Edit->Material Coordinates…, if the model shown is not MatCoordStandard, hit ‘Cancel’ in that window. Then within the Continuity client, go to File→Library→Search… and enter MatCoordStandard in the top search window. Right-click on the model that was found and select Load. It is IMPORTANT that you select the third option in the pop-up windows that starts with Retain current problem…. Then, if you go back to Mesh->Edit->Material Coordinates, you should now have the Standard Material Coordinate model equations and parameters. Compile the model under Compile…→Compile Code as: C Double Precision. When the compilation has completed, check the Parameters to make sure that the fiber_angle parameter is set to the Fiber Angle field. You can change the values of the parameters later without having to recompile the model. Click OK to exit the model editor.

• Continuity Command Summary Mesh→Edit→Material Coordinates… Confirm that MatCoordStandard model is listed If MatCoordStandard is not listed:    File→Library→Search…    Search for MatCoordStandard    Load model and select third option in Resolve Conflict form Mesh->Edit->Material Coordinates… Compile…→Compile Code as: C Double Precision Check the Parameters to make sure that the fiber_angle parameter is set to the Fiber Angle field Click OK Save model

   

Solving a Simple Biomechanics Problem in Prolate Spheroidal Coordinates

Mesh Refinement

To solve a Biomechanics problem, we will first need to refine the mesh to allow for reasonably converged stress and strain solutions. Here we will refine the mesh into 5 elements longitudinally and 2 transmurally. Select Mesh→Refine… and under “New Element per old element in:” select 1 for “xi1”, 5 for “xi2” and 2 for “xi3”. Click OK. Then re-render the mesh after first remembering to use the Mesh→Calculate Mesh… command. Your mesh should now have ten elements. Note that the node numbers will have changed but the labels will be preserved. More importantly the new nodes will have inherited the labels in a logical way. The new node at the apex will have the label APEX, the new nodes at the base will have the label BASE, the new nodes along the inside surface will be labeled ENDO and the new nodes along the outside will be labeled EPI. Use the Edit→Nodes… command from the Mesh menu to open the Nodes form and confirm this.

• Continuity Command Summary Mesh→Refine… Select 1 for “xi1”, 5 for “xi2” and 2 for “xi3” under “New Element per old element in:” Click OK Calculate Mesh Render Elements Save model

   

Biomechanics Problems

The Biomechanics module of Continuity uses a Newton iterative method to solve the continuum equations of motion for nonlinear, anisotropic, large deformation elasticity problems subject to traction and/or displacement boundary conditions. You can load the Biomechanics menu and module using the add module toolbar button, or from File->Load Modules. Various other special features of biomechanics problems, notably viscoelasticity, dynamic active force generation and boundary conditions generated by circulatory models can be included too. The necessary governing equations include the nonlinear kinematic strain-displacement relations for large deformations, conservation of linear and angular momentum equations from continuum mechanics, a nonlinear constitutive equation for hyperelastic materials, and in the case of incompressible models, the Lagrangian mass conservation law. The weighted residual formulation of the finite element equations uses the Galerkin finite element method, which is equivalent in this case to solving the virtual work integral equations of mechanics.

We will also change our basis functions for the undeformed and deformed geometric variables from tri-linear to linear-bicubic to allow a more continuous and higher order variation of stress and strain transmurally and longitudinally. The refine command will automatically generate the correct derivatives we need for the undeformed geometry. If you still have the Nodes form open, change the basis functions for “Coordinate 1”, “Coordinate 2” and “Coordinate 3” to Linear-Cubic-Cubic Hermite 3*3*3. Click OK. Observe that the non-zero parameters that the mesh refinement has calculated and inserted into the form coincide to the derivatives of the coordinates with respect to the directions in which they are changing, i.e. Coordinate 1 (Lambda) is changing in the xi3 direction, Coordinate 2 (mu) is changing wrt xi2.

• Continuity Command Summary Mesh→Edit→Basis… and select Hermite Basis Function→3D→Linear-Cubic-Cubic Use 3 integration/collocation points for each local Xi coordinate. Click Add then OK to submit. Mesh→Edit→Nodes Change Coordinates 1, 2, and 3 to be Linear-Cubic-Cubic Hermite 3*3*3 Click OK

   

Boundary Conditions

We will use the Biomechanics Edit→Boundary Conditions menu choice to assign boundary conditions. But first use the Biomechanics→Update→Initial conditions from undeformed nodes command to copy the undeformed nodal coordinates and basis functions to the initial conditions tab of the Boundary Conditions form. Next, we will use the Edit→Boundary Conditions form to assign nodal displacement boundary constraints. Click on the Deformed Coordinate 1 tab to specify constraints on the deformed Lambda coordinate, the Deformed Coordinate 2 tab to specify constraints on the deformed Mu coordinate, and the Deformed Coordinate 3 tab to specify constraints on the deformed Theta coordinate. For the node numbers, you can use the labels BASE.*, APEX.* and BASE_EPI to specify the list of nodes for each constraint as shown in the table. (note the “.*” is a wildcard option that will capture the nodes labeled BASE_EPI, APEX_ENDO etc).

Next we can assign pressure boundary conditions using the External Pressure tab of the Boundary Conditions Form. In the “Element Number:” field enter a comma-separated list of the element numbers on the inner surface of the model. These are elements 1-5. (You can use the tools in View→Show→OpenMesh to confirm this). Next to “Select Pressure Type”, select “Incremental” and assign the pressure increment to the inner faces. For passive filling of the left ventricle 0.1 kPa should work well.

• Continuity Command Summary Edit→Load Modules→Biomechanics Biomechanics→Update and select Initial conditions from undeformed nodes Biomechanics→Edit→Boundary Conditions Click on Deformed Coordinate 1 to specify constraints on the deformed Lambda coordinate Click on Deformed Coordinate 2 to specify constraints on the deformed Mu coordinate Click on Deformed Coordinate 3 to specify constraints on the deformed Theta coordinate Click on External Pressure and define inner surface elements Select “Incremental” under “Select Pressure Type” Specify pressure for inner faces Click OK

   

Constitutive Model

To specify the material properties, we use the Edit→Constitutive Model command from the Biomechanics menu to open the Model Editor. Here you can edit an existing model or author and compile a new one provided the compilation environment is set up. One model that is in the Continuity database (BM_Ortho_of_Lagrangian_strains_comp_U8_sympy) is orthotropic exponential slightly incompressible strain energy function of Lagrangian strains proposed by Usyk and McCulloch for canine myocardium in which the strain energy is a function of strain components referred to local fiber-sheet-normal material axes. Of course, without recompiling we can change the parameters, of which there are seven, an overall scaling coefficient and one coefficient for each of the six strain components in fiber-sheet-normal coordinates (E_ff, E_ss, E_nn, E_fs, E_fn, E_sn) that appears in the exponent of this strain-energy function. However, we will start with a slightly simpler transversely isotropic exponential strain energy function in which the fiber direction of the tissue is stiffer than the transverse direction but there is no difference between the sheet and sheet-normal (cross-fiber) directions. This can be found in the Continuity database by searching for: BM_TI_of_Lagrangian_strains_comp_sympy), and can be loaded by right-clicking and selecting load and the selecting the third option in the Resolve Conflict window. This model has only five parameters instead of eight for the orthotropic model. The last parameter is a “Bulk Modulus” of the “slightly incompressible” term in the model. It is relatively high and penalizes changes in volume as the tissue pressure rises. The higher this number, the less the volume of the tissue will change as the material deforms. A number that is too low can make the solution sensitive to the specific value. As the value gets higher the solution changes less but can be harder to obtain. We will start with the default transversely isotropic exponential model. Later we will switch to the orthotropic model. Leave the parameters in their default values. Click Submit.

• Continuity Command Summary File→Library→Search Search for “BM_TI_of_Lagrangian_strains_compy_sympy” Right-click on found model and select “Load” From Resolve Conflict window, select third option (“Retain current problem…”) Biomechanics→Edit→Constitutive Model… Verify model is loaded and click “OK” to exit

   

Check Model Variables Size Allocation

Before we solve the biomechanics problem, we should check that the dimensions of all the array variables in the model are correct to avoid simple computational errors, particularly indexing errors. Go to View→Edit Dimensions and ensure that All Dimensions are OK is displayed. If there is a message about variable sizes being too small or too large, click Apply Marked Recommendations to automatically change the allocated sizes of all incorrectly sized variables. If you click Show Details, you can view each variable and its current and recommended sizes.

• Continuity Command Summary View→Edit Dimensions Verify all dimensions or OK, or if not, accept the recommended corrections by selecting “Apply Marked Recommendations” Click “OK” to exit Do a send Calculate Mesh Note:, if you do a send and get an error about surface_origin. Save your model, then restart Continuity and reload your model. If you do send after reloading your model, it will work fine.

   

Solving the Biomechanics Model and Displaying the Results

Now all that remains is to issue the Biomechanics→Solve command. If you specify ten time steps of 1.0 units each, then 10 previously specified pressure increments will be applied; the final pressure will be 1.0 kPa. Note that the solve form includes stopping criteria for the iterative solver based on the sum of the solution vector increments, the sum of the unconstrained residuals and the maximum number of iterations. For this problem, the default setting will work, but that is not always the case. If the solution is having difficulty converging, you can relax the stopping criteria by increasing the number of iterations and the tolerance for the sum of solution increments and sum of unconstrained residuals (usually up to 1.00e-002). Decreasing the time step by orders of magnitude is also helpful.

• Continuity Command Summary Biomechanics→Solve Nonlinear Change “Number of Steps” to 10 Click “Solve”

   

Note that the next command makes use of a Field Variable which must have a basis function assigned in the Nodes Form. For Field Variable 1 in the Field Vector 1 tab of the Edit→Nodes form, choose Linear-Linear-Linear Lagrange 3*3*3. Once the problem has solved we can use the Biomechanics→Render→Surface command: under “Variables” select T- Cauchy Stress Tensor to visualize individual components of the stress or strain tensor solutions on specified element surfaces. If you want to change the color, go to View→Show Open Mesh …, Select Textured Field from the Rendering order list, then Click on Colors tab select RGB or BGR for Color Scheme:, then Close the window. If you specify ξ₁₁=0 or ξ₁₁=0.5 you will be able to see results rendered on longitudinal planes. Note that fiber strain (xx) is more uniform than the transverse (yy or zz) components.

• Continuity Command Summary Mesh→Edit→Nodes On the Fiber Angles tab, make sure that Linear-Linear-Linear Lagrange 3*3*3 is listed under Fiber Angle Click “OK” Biomechanics→Render→Surface select variable T- Cauchy Stress Tensor Click “OK”

   

Now save your model file, restart Continuity, reload the model and change the Constitutive model to the orthotropic exponential model. Change the Bulk Modulus in the parameter tab from 1000 to 400. For this model to converge, you will need to change some of the solve settings. Change the time step to 0.5 and take 20 steps instead of ten. Reduce Delta in the solution scheme tab from 1.e-4 to 1.e-7. Compare the solutions with the previous case.

Prebuilt models

Prebuilt models ready for the biomechanics solve can be found in the Continuity model database. To access them, go to File→Library→Search… and find GettingStartedBMtrans.cont6 for the transversely isotropic case and GettingStartedBMortho.cont6 for the orthotropic case. To load them, simply right-click on the title of the model and select Load.

Fitting

• Getting started with fitting

Electrophysiology

• Getting started with Electrophysiology

• Getting started with Electrophysiology on a cluster

• Getting started with Electrophysiology (MFHN model)

Biomechanics

• Getting started with Biomechanics

Slide Show

See a slide show with all the images from this article.