# Getting Started With Continuity

Introductory Slides YouTube Video Tutorial

Contents

- Getting Started With Continuity
- Slide Show

## 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_{1} and u_{2}) of the function at the nodes, and thus the basis functions (φ_{1} and φ_{2})are the corresponding linear functions that when multiplied by u_{1} and u_{2}, respectively, and added together result in the linear interpolation between u_{1} and u_{2}. 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_{0}) 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 ξ_{1} and ξ_{2} and 3-D hexahedral trilinear Lagrange elements as functions of ξ_{1}, ξ_{2}, and ξ_{3}.

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. ξ_{1}), "Xi 2" (i.e. ξ_{2}), and "Xi 3" (i.e. ξ_{3}). 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 (ξ_{1},ξ_{2},ξ_{3})=(0,0,0) and ending with the node at (ξ_{1},ξ_{2},ξ_{3})=(1,1,1). By convention for meshes defined in curvilinear coordinates, we make ξ_{1} coincide (to the extent possible) with the polar (θ) axis, ξ_{2} with the "longitudinal" (z, φ or μ) axis, and ξ_{3} 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 ξ_{1} 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, ξ_{1} 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 (ξ_{1}), then longitudinally (ξ_{2}) from apex to base, and finally radially (ξ_{3}) 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 *h*x*a*, which is 0.45x4.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_{1}) 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_{1},Y_{2},Y_{3}); curvilinear World coordinates (Θ_{1},Θ_{2},Θ_{31}); local finite element coordinates (ξ_{1},ξ_{2},ξ_{3}); and local *material* coordinates (x_{1},x_{2},x_{3}). 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 ξ_{1}, one in the ξ_{1}-ξ_{2} 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:

**MatCoordStandard**

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 ξ_{11}=0 or ξ_{11}=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

## Electrophysiology

## Biomechanics

# Slide Show

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