Material models

Material models#

When we solve the mechanical model, we are solving for an equilibrium of stresses. Here, both active and passive tension plays in, as well as boundary conditions. The passive tension is determined by a strain energy function, which again depends on certain material parameters.

In an earlier demo we demonstrated how to perform virtual stretch and shear experiments. In this repo, we’ll redo the stretch experiment – and compare two material models; one Guccione kind of model, and one Holzapfel kind of model. Note that the Holzapfel model is fitted with experimental data, while the Guccione model is simply adjusted to be “about as stiff”. Both assume a transversely isotropic formulation in the intracellular subdomain, and an isotropic formulation for the extracellular subdomain.

We will use the code from the main repository emimechanicalmodel, which is based on FEniCS. We can import all libraries we need in this demo, and load the meshes and subdomains:

import numpy as np
import matplotlib.pyplot as plt

import fenics as f
import emimechanicalmodel as emi_m

f.set_log_level(30)        # less verbatim output from the Newton solver compared to default

mesh_file = "tile_connected_10p0.h5"
mesh, volumes = emi_m.load_mesh(mesh_file)

Mesh and subdomains loaded successfully.
Number of nodes: 335, number of elements: 1364

The strain energy functions#

When we solve the mechanical problem, we want to find a displacement field \(\mathbf{u}\) which gives us an equilibrium of stresses. We can write

\[\Psi = \Psi_{active} + \Psi_{passive} + \Psi_{incompressible}\]

and then define \(\mathbf{P}\) based on \(\Psi\) and the deformation tensor \(\mathbf{F}\). See the previous chapter (or the papers) for the mathematical details.

In this chapter we’ll focus on \(\Psi_{passive}\), which we can separate it into two separate contributions:

\[\begin{split} \Psi_{passive} (F) = \begin{cases} \Psi_{i} (F) \qquad x \in \Omega_i \\ \Psi_{e} (F) \qquad x \in \Omega_e \end{cases} \end{split}\]

Guccione model#

We have adapted a traditional Guccione model (or a Fung model) to the EMI framework, by defining two strain energy functions, one for each subdomain. For this, we use components of the Cauchy-Green strain tensor \(\mathbf{E}\).

\[ \Psi_i = C_i (e^{Q_i} - 1)\]

The parameter \(C_i\) has unit kPa, and adjusting will adjust the stiffness of the cell. Furthermore, the term \(Q_i\) determines the exponential behavior of increasing amounts of strain:

\[ Q_i = b_{f, i} E_{11}^{2} + b_{t, i}(E_{22}^{2} + E_{33}^{2} + E_{23}^{2} + E_{32}^{2}) + b_{ft, i}(E_{12}^{2} + E_{21}^{2} + E_{13}^{2} + E_{31}^{2}) \]

The exponential (“b”) parameters determine the ratio between the fiber direction (\(b_{f, i}\)), the transverse direction (\(b_{t, i}\)) and the transverse direction (\(b_{ft, i}\)). If they are all equal, the material can be called isotropic.

Next, for the extracellular space, we will use

\[ \Psi_e = C_e (e^{Q_e} - 1)\]

where

\[ Q_e = b_e (E_{11}^{2} + E_{22}^{2} + E_{33}^{2} + E_{23}^{2} + E_{32}^{2}) + E_{12}^{2} + E_{21}^{2} + E_{13}^{2} + E_{31}^{2}) \]

Mathematically, this term is equal to \(Q_i\) but with all parameters equal to \(b_e\). Mechanically, this implies that the material is isotropic with material properties determined by \(C_e\) and \(b_e\).

Holzapfel model#

As for the Guccione kind of model, we define separate strain energy functions \(\Psi_i\) and \(\Psi_e\). The Holzapfel model is invariant-based. It is, arguably, less advanced than a traditional Holzapfel model, in terms of only taking into account two.

We first define the right Cauchy-Green tensor (in fact, \(\mathbf{E}\) is usually defined based on \(\mathbf{C}\))

\[ \mathbf{C} = \mathbf{F}^T \mathbf{F} \]

then we have

\[\begin{split} I_1 = \mathrm{tr} (\mathbf{C}) \label{eq:intro:I_1} \\ \end{split}\]

and

\[\begin{split} I_{4, f} = \mathbf{f} \cdot (\mathbf{C} \mathbf{f}) \label{eq:intro:I_4e1} \\ \end{split}\]

where \(I_{4, f}\) simply quantiifes the change in the fiber direction. Next, we define

\[ \psi(\mathbf{F}) = \frac{a_i}{2 b_i} (e^{b_i (I_1 - 3)} - 1) + \frac{a_{if}}{2 b_{if}} (e^{b_{if} \| I_{4f} - 1 \|_{+}^2} - 1) \]

where \(\| \cdot \|_{+}\) is a conditional term given by \(\| I \|_{+} = \mathrm{max}(I, 0)\) (i.e., 0 if \(I\) is below zero). For the extracellular subdomain, we define

\[ \psi_{e}(\mathbf{F}) = \frac{a_e}{2 b_e} (e^{b_e (I_1 - 3)} - 1). \]

Stretching and active contraction#

We have two material models – both based on well-known, widely used traditional tissue-level models. They are fairly interchangable, but not equal – note that this is, in parts, due to that the Guccione model has not been parametrized in a rigorous way. We’ll compare stretch in different directions and contraction values, tracking load values for the stretching experiments and intracellular fiber direction strain and stress values for the contraction experiments.

material_models = ["guccione", "holzapfel"]
experiments = ["stretch_ff", "stretch_ss", "stretch_nn"]

stretch_values = np.linspace(0, 0.15, 15)

fig, axes = plt.subplots(1, 3, figsize=(15, 3))

for material_model in material_models:

    for experiment, axis in zip(experiments, axes):
        emimodel = emi_m.EMIModel(
                mesh,
                volumes,
                experiment=experiment,
                material_model=material_model,
            )

        load_values = np.zeros_like(stretch_values)

        # then run the simulation
        for step, s_str in enumerate(stretch_values):
            emimodel.assign_stretch(s_str)
            emimodel.solve()

            load_values[step] = \
                emimodel.evaluate_normal_load()

        axis.plot(100*stretch_values, load_values)

axes[0].set_xlabel("Stretch (%)")
axes[1].set_xlabel("Stretch (%)")
axes[2].set_xlabel("Stretch (%)")
axes[0].set_ylabel("Load (kPa)")

axes[1].legend(["Guccione", "Holzapfel"])

plt.show()

Calling FFC just-in-time (JIT) compiler, this may take some time.

Domain length=102.0, width=20.0, height=24.0

0.0
Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

load:  0.0
0.010714285714285714

load:  0.9647807287871485
0.02142857142857143

load:  1.9977946039552554
0.03214285714285714

load:  3.128628313776314
0.04285714285714286

load:  4.392473145333273
0.053571428571428575

load:  5.832412213998076
0.06428571428571428

load:  7.502369161843769
0.075

load:  9.470986969463286
0.08571428571428572

load:  11.826865874387954
0.09642857142857143

load:  14.685770709977486
0.10714285714285715

load:  18.200670548862536
0.11785714285714285

load:  22.575858738118253
0.12857142857142856

load:  28.086992615592067
0.1392857142857143

load:  35.10979744251787
0.15

load:  44.16157129429943
Domain length=102.0, width=20.0, height=24.0

0.0
Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

load:  0.0
0.010714285714285714

load:  0.3861671891962323
0.02142857142857143

load:  0.779329375028315
0.03214285714285714

load:  1.183722091936638
0.04285714285714286

load:  1.6036606504367696
0.053571428571428575

load:  2.043627400655194
0.06428571428571428

load:  2.5083692244072
0.075

load:  3.0030053857557455
0.08571428571428572

load:  3.533146221302221
0.09642857142857143

load:  4.105024253222878
0.10714285714285715

load:  4.725640686741672
0.11785714285714285

load:  5.402931633869478
0.12857142857142856

load:  6.145959727330364
0.1392857142857143

load:  6.965138130302515
0.15

load:  7.872495445631977
Domain length=102.0, width=20.0, height=24.0

0.0
Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

load:  0.0
0.010714285714285714

load:  0.3202565185579636
0.02142857142857143

load:  0.6450736575890885
0.03214285714285714

load:  0.9769105571831688
0.04285714285714286

load:  1.3182440249821743
0.053571428571428575

load:  1.671600179951388
0.06428571428571428

load:  2.039589432514045
0.075

load:  2.42494507242055
0.08571428571428572

load:  2.830565551029815
0.09642857142857143

load:  3.259560545878559
0.10714285714285715

load:  3.7153010452269037
0.11785714285714285

load:  4.201473934456764
0.12857142857142856

load:  4.7221418566678155
0.1392857142857143

load:  5.281809425907178
0.15

load:  5.885497181095201
Domain length=102.0, width=20.0, height=24.0

0.0
Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

load:  3.286738961392422e-16
0.010714285714285714

load:  0.4134900403593821
0.02142857142857143

load:  0.8844810028495402
0.03214285714285714

load:  1.4439205180459815
0.04285714285714286

load:  2.131952291123514
0.053571428571428575

load:  3.0055951255903213
0.06428571428571428

load:  4.1475331179450805
0.075

load:  5.680252004114171
0.08571428571428572

load:  7.789105904426988
0.09642857142857143

load:  10.760421435736967
0.10714285714285715

load:  15.045687831721315
0.11785714285714285

load:  21.371894821558612
0.12857142857142856

load:  30.9350278389588
0.1392857142857143

load:  45.746416239068026
0.15

load:  69.26606752990033
Domain length=102.0, width=20.0, height=24.0

0.0
Calling FFC just-in-time (JIT) compiler, this may take some time.

load:  2.326665247791059e-16
0.010714285714285714

load:  0.09897749001031107
0.02142857142857143

load:  0.20292674531579205
0.03214285714285714

load:  0.3157854901762425
0.04285714285714286

load:  0.4406997969488145
0.053571428571428575

load:  0.5806949589561573
0.06428571428571428

load:  0.7392142182674761
0.075

load:  0.9204508482058973
0.08571428571428572

load:  1.1295974471162054
0.09642857142857143

load:  1.3731034222538472
0.10714285714285715

load:  1.658990109070843
0.11785714285714285

load:  1.9972564491679305
0.12857142857142856

load:  2.400406843250593
0.1392857142857143

load:  2.8841392880146515
0.15

load:  3.468243665395189
Domain length=102.0, width=20.0, height=24.0

0.0
Calling FFC just-in-time (JIT) compiler, this may take some time.

load:  2.22520192647767e-16
0.010714285714285714

load:  0.07374225734707336
0.02142857142857143

load:  0.15060874189786755
0.03214285714285714

load:  0.23316350928522817
0.04285714285714286

load:  0.3236146409226338
0.053571428571428575

load:  0.4238985299860558
0.06428571428571428
load:  0.535862287384527
0.075

load:  0.6614630109977988
0.08571428571428572

load:  0.8029433802286245
0.09642857142857143

load:  0.962978897505429
0.10714285714285715
load:  1.1448092435546937
0.11785714285714285

load:  1.3523701180197583
0.12857142857142856

load:  1.5904405133806538
0.1392857142857143

load:  1.8648183207026894
0.15

load:  2.182536257715045

../_images/3bcb32e982e732e865e744729568e54f6209a6a3b34e7656cd59b21631fc0699.png

Contraction#

time = np.linspace(0, 500, 125)       # 500 ms with 125 steps
active_precomputed = emi_m.compute_active_component(time)

material_models = ["guccione", "holzapfel"]

fig, axes = plt.subplots(1, 2, figsize=(15, 3))

for material_model in material_models:
    emimodel = emi_m.EMIModel(
            mesh,
            volumes,
            experiment="contraction",
            material_model=material_model,
        )

    fiber_dir_strain_i = np.zeros_like(active_precomputed)
    fiber_dir_stress_i = np.zeros_like(active_precomputed) 
    subdomain_id = 1         # intracellular subdomain

    # then run the simulation
    for step, a_str in enumerate(active_precomputed):
        emimodel.update_active_fn(a_str)
        emimodel.solve()

        fiber_dir_strain_i[step] = \
            emimodel.evaluate_subdomain_strain_fibre_dir(subdomain_id)
        fiber_dir_stress_i[step] = \
            emimodel.evaluate_subdomain_stress_fibre_dir(subdomain_id)

    axes[0].plot(time, fiber_dir_strain_i)
    axes[1].plot(time, fiber_dir_stress_i)

axes[0].set_xlabel("Time (ms)")
axes[1].set_xlabel("Time (ms)")
axes[0].set_ylabel("Strain (-)")
axes[1].set_ylabel("Stress (kPa)")

axes[0].set_title(r"Fiber direction strain, over $\Omega_i$")
axes[1].set_title(r"Fiber direction stress, over $\Omega_i$")

axes[1].legend(["Guccione", "Holzapfel"])

plt.show()

Domain length=102.0, width=20.0, height=24.0

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Calling FFC just-in-time (JIT) compiler, this may take some time.

Domain length=102.0, width=20.0, height=24.0

../_images/fb28a9c730b4d35dd161ddf450dafa0ae4dd7f2492ac4a78bd5d6910acda2389.png