Remodeling-induced stiffening of fibrous networks under deformation

University of Amsterdam

MSc Physics

Physics of Life and Health

Master Thesis

Remodeling-induced stiffening of

fibrous networks under

deformation

by

Luka Liebrand

5898617

October 2015 Credited 60 ECTs Research conducted between: December 2014 - October 2015

Supervisors:

Prof. dr. Fred MacKintosh

Dr. Mahsa Vahabi

Examiners:

Prof. dr. Fred MacKintosh

Dr. Edan Lerner

Abstract

The mechanical properties of disordered extracellular biopolymer networks play an important role in the elasticity and stability of tissue. Studies have shown that these networks exhibit nonlinear stiffening when subjected to strain. There is a linear stiffness regime for small shear deformations, while upon reaching higher levels of deformation these networks exhibit dramatic nonlinear stiffening of at least an order of magnitude. This reversible strain-stiffening can already be explained by simple athermal network models of fibers with linear stretching and bending interactions.

Recently, it has been found that fibrin networks exhibit additional stiffening behavior that is non-reversible. Networks of fibrin, which is a biopolymer that is necessary for blood clotting, are stiffer after being compressed. Optical tweezer experiments also show that fibrin fibers are capable of forming new bonds when brought into contact with each other. In this study, it is investigated whether remodeling, in the form of forming additional bonds within the network, can lead to network stiffening by developing a model with this type of remodeling.

We present an athermal fibrous network model that allows forming new bonds. Simulations wherein networks were compressed and decompressed show that remodeling leads to an increased stiffness, with the same qualitative be-havior as observed in fibrin. Additionally, a previously unknown power-law relationship was found between the residual normal stress and the increase in the shear modulus after compression. For both experiments and simulations, the obtained exponent was approximately two thirds.

Nederlandse samenvatting

De meeste levende systemen zijn complex en hebben zich onder druk van natu-urlijke selectie gespecialiseerd in het uitvoeren van een specifieke functie. Zo hebben natuurlijke weefsels, afhankelijk van hun rol in het lichaam, verschil-lende interessante mechanische eigenschappen. Deze mechanische eigenschap-pen worden in grote mate ontleend aan complexe netwerken van biologische eiwitten. Een goed voorbeeld van een complex biologisch weefsel met specifiek ontwikkelde eigenschappen is de huid, die voornamelijk bestaat uit collageen. De huid is elastisch en geeft makkelijk mee bij kleine vervormingen, maar wordt opeens zeer stug bij grote mate van deformatie. Deze eigenschap houdt de huid enerzijds soepel en zorgt ervoor dat we makkelijk kunnen bewegen, maar zorgt er anderzijds voor dat de huid niet te vervormd raakt en het onderliggende weefsel kan blijven beschermen.

Dat biologische weefsels stijver zijn op het moment dat ze worden vervormd is algemeen bekend. De natuurkundige oorsprong van deze omkeerbare versti-jving ligt in de eerder genoemde netwerken van biologische eiwitten, maar de kennis van deze netwerken is nog beperkt. Op dit moment zijn er onderzoeken gaande die ervoor moeten zorgen dat we de mechanische eigenschappen van eiwitnetwerken en de onderliggende fysica beter te begrijpen. Een interessant eiwit is fibrine. Fibrine is medeverantwoordelijk voor het stollen van bloed, do-ordat het een nauw netwerk vormt dat bloedplaatjes vangt en de bloedstroom stopt. Recent onderzoek heeft namelijk aangetoond dat netwerken van fibrine niet alleen stijver zijn op het moment dat ze zijn opgerekt, maar ook als ze eerst zijn samengeperst en daarna weer tot hun originele grootte zijn teruggebracht. De onomkeerbare verstijving van fibrinenetwerken heeft dit onderzoek gemo-tiveerd. Met behulp van een netwerkmodel wordt getracht deze verharding uit te leggen. Omdat er tot dusver geen netwerkmodel dat dit fenomeen kan verklaren beschikbaar was, werd er voor dit onderzoek een ontwikkeld. Het mechanisme dat hierin verantwoordelijk is voor het verstijven komt overeen met een an-dere recente observatie van fibrine. Individuele fibrinevezels kunnen namelijk nieuwe verbindingen vormen als ze dicht genoeg bij elkaar worden gebracht. Dit plakkende effect zorgt ervoor dat er minder bewegingsvrijheid is voor de vezels in het netwerk, waardoor het netwerk uiteindelijk stijver wordt.

Het gepresenteerde model bestaat uit een onregelmatige ori¨entatie van lineair elastische componenten. Dit wil zeggen dat de kracht proportioneel is met de uitrekking of buiging van elk onderdeel. Het is al aangetoond dat een combinatie van lineaire componenten op netwerkniveau tot non-lineair gedrag kan leiden. Om de onomkeerbare verstijving aan het netwerkgedrag toe te voegen werd het vezels toegestaan om aan elkaar te binden als ze binnen een bepaalde afstand kwamen. Tijdens de compressie die op het netwerkmodel werd uitgeoefend in de simulaties, kwamen de vezels binnen deze bindingsafstand en werden ze dus verbonden. Het resultaat was een verstijving van het netwerkmodel, met kwalitatief gelijk gedrag als fibrine. Daarnaast kan, door de bindingsafstand te vari¨eren, de kracht van deze verstijving be¨ınvloed worden. We hebben ook een relatie gevonden die de uiteindelijke verstijving van het netwerk koppelt aan de druk die nodig is om het netwerk opgerekt te houden. Het natuurkundig principe dat hieraan ten grondslag ligt is echter nog niet bekend en vereist verder onderzoek.

Contents

1 Introduction 4

2 Theory 8

2.1 Microscopic interactions . . . 8

2.2 Macroscopic interactions: Connected fibers . . . 10

2.3 Choosing the network structure . . . 11

2.3.1 Maxwell’s criterion . . . 11

2.3.2 Creating the network from a lattice structure . . . 12

2.3.3 Midpoints . . . 13

2.3.4 From a triangular lattice to a disordered network . . . 13

2.3.5 Boundary conditions . . . 14

2.4 Extracting mechanical properties . . . 15

2.4.1 Measuring mechanical properties . . . 15

2.4.2 The rheometer . . . 16

2.4.3 Energy and deformations in simulated networks . . . 16

2.4.4 Mechanical properties in simulated networks . . . 19

2.5 Behavior of fibrous networks . . . 20

2.5.1 Strain-stiffening in biopolymer networks . . . 20

2.5.2 Fibrin network stiffening . . . 21

2.6 Does remodeling explain fibrin network stiffening? . . . 23

2.6.1 Probing fibrin fiber stickiness with optical tweezers . . . . 23

2.6.2 Adding remodeling to the network simulations . . . 24

3 Results 29 3.1 Network stiffening as a result of cyclic deformation . . . 29

3.1.1 The effect of the remodeling distance on network stiffening 31 3.2 A power law relates network stiffness to the normal stress . . . . 32

4 Discussion and Conclusion 36

Chapter 1

Introduction

Living systems are often complex in their structure and dynamics. In order to fully understand how such a living system functions, one has to completely understand their constituent parts and their interactive behavior. Therefore, it can be beneficial to isolate the individual components of the system of interest and study them separately. These components range from one type of tissue to complete organs. While the workings of a wide variety of living systems are already quite well known, others remain elusive; one of the open questions is how biopolymers and their networks relate to tissue characteristics.

Many soft tissues in living organisms derive their mechanical properties from disordered networks of different biopolymers [1, 2]. These protein networks do not only passively give structure to their surroundings, but may also actively play a role in their dynamics, as in the case of the acto-myosin network in-side cells [3, 4]. In figure 1.1, two examples of biopolymer networks inin-side the body are shown: acto-myosin and collagen. Acto-myosin is responsible for force generation within cells and controls a cell’s shape, whereas collagen connects elements in the extra-cellular matrix and promotes inter-cellular stability.

While the benefits of the elasticity found in tissue are quite clear, the physics that govern the mechanical properties of disordered biopolymer networks are not yet well-understood. That is why these are currently studied both ex-perimentally and theoretically. Because of the complexity of biological tissue, experimental studies often apply a bottom-up approach to recreate networks, by reconstituting networks from an isolated type of polymer [4, 5, 6, 7]. This practice also allows highlighting differences between different biopolymer types, where variations in mechanical properties can be linked to their function inside the body.

Complementary to experiments, work is being done on theoretical models whose aim is to explain or predict the mechanical properties of biopolymer net-works [1, 8], to eventually gain understanding of the underlying physical prin-ciples. In these models, how the individual fibers are connected into a network structure is as important as the properties of the individual fibers. Arguably the most important characterizing parameter of the network structure is the average connectivity of the network, z. This parameter indicates the average coordination number (i.e. the number of other nodes one node is connected to) of nodes in the network. A network with a higher z has fewer degrees of freedom (DOF) and that is why it is related to the network’s stability and resistance to

Figure 1.1: Confocal fluorescence micrographs of (a) a fibroblast (yellow) in an extracellular collagen network structure (green), and (b) an acto-myosin network, with nuclei (cyan), microtubules (yellow) and actin filaments (pur-ple). Images from: http://biochem.wustl.edu/confocor/collagenbar.png and http://publications.nigms.nih.gov/insidelifescience/images/microtubule big.jpg, respectively.

Figure 1.2: The nonlinear stiffness K of collagen under shear strain γ is shown for measurements with different concentrations (data points) along with a the-oretical prediction (dashed lines). The inset shows the relationship between the concentration and the dimensionless bending rigidity ˜κ (see section 2.1). Figure from Sharma et al. [16].

Figure 1.3: Shear modulus (G) vs axial compression () for a reconstituted fibrin gel under (cyclic) axial compression. Each successive compression cycle reaches a higher , resulting in a higher G each time the sample is decompressed to = 0. Data (unpublished) by: Bart Vos et al.

deformation. This was already understood by Maxwell [9]: According to his theory, a network of Hookean springs is only stable if the average connectivity is at least two times the dimensionality of the network (z ≥ 2d). Because the connectivity in tissue is far below Maxwell’s connectivity treshold of z = 6, ef-fects beside stretching must also play a role in network stability. For example, z≈ 3.4 in collagen networks [10].

Earlier work has shown that there are multiple ways to account for the stability of sub-marginal (where z< 2d) networks [11, 12, 13]. For example, ex-ternally applied strain can be responsible for stability in sub-marginal networks, as well as bending interactions between fiber segments [11]. In fact, a network of Hookean springs with bending interactions can already exhibit similar mechan-ical properties to real biopolymer networks, such as nonlinear strain-stiffening [14, 7, 15]. A comparison between a theoretical prediction and experimental values of the shear modulus of collagen networks under shear strain by Sharma et al. [16] is shown as an example in figure 1.2. Here, the shear modulus, which is defined as the ratio of the shear stress to strain, indicates the stiffness of the network. Clearly visible is the strain-dependent nonlinear stiffening; this stiffening is present in all biological tissues, and can be explained using thermal [7] and athermal models [16, 11]. In athermal models, the source of the non-linear stiffening is purely mechanical, as the fibers are less resistant to bending than stretching, so bending is energetically favorable. Hence, there is a more compliant regime at small deformations, where the network can cope by merely bending, and a stiffer regime when the fibers are stretched.

this field, with (rheometer) experiments motivating theoretical research and vice-versa [4]. Actually, this study is motivated by recent rheometer experi-ments done on fibrin, a biopolymer that is important in blood clotting. These experiments have shown that fibrin networks show irreversible stiffening after being subjected to compression. This is expressed in figure 1.3, where a com-pressed and decomcom-pressed fibrin sample shows a much larger shear modulus than initially.

Since the current models have not been developed to explain fibrin’s irre-versible stiffening behavior, a new model is needed to account for this effect. One possible explanation for the stiffening is that during deformation fibrin networks form new bonds, which act like molecular motors. Molecular motors can generate forces to locally strain the network. In this way, a fibrin network ends up stressed, even when it is eventually decompressed. Such a remodeled network has a lower number of DOF and a higher resistance to deformation.

Evidence that fibrin fibers can stick together and form new bonds has re-cently been found with optical tweezer experiments that are presented in section 2.6.1. However, proof that the observed stiffening can be attributed to bond formation is still required. This study seeks to answer the question whether bond formation can lead to stiffening in fibrous networks, and to what extent. This question will be answered with the help of a newly-developed model that allows for remodeling during deformation. Finally, the results will be compared with those of fibrin networks.

Chapter 2

Theory

In order to understand all aspects of a network, one must first look at its con-stituent parts. The networks that are considered in this thesis are made up of fiber segments connected at junctions where they can hinge called nodes. As mentioned, these fiber segments have elastic properties, and a finite resistance to bending. In the following sections the microscopic interactions (i.e. interac-tions on the segment-level) within the network will be explained, as well as how they lead to the macroscopic behavior of the network as a whole.

2.1

Microscopic interactions

Biopolymer networks can have different elastic properties that depend on fiber types and network structure. The main reason for the fiber’s elasticity is to prevent structural damage to the tissue upon deformation [17, 18]. The elas-ticity of single fibers is linear for small deformations, but becomes non-linear for increasing strain levels [19, 20], thus ensuring that the tissue does not get displaced too far.

The linear elastic behavior of a biopolymer fiber can be captured by a Hookean spring. Because it is possible to understand the nonlinear behavior of biopolymer networks with only linear components [11], Hookean springs rep-resent fiber segments in this work. The Hamiltonian of a Hookean spring is, of course, well-known. Throughout this thesis the following form is used:

Hstretch= 1 2 µ l0(δl) 2_{, (2.1)}

where µ is the elastic modulus and δl is the extension with respect to the rest length of the spring, l0. As can be seen from this formula, the distance between

the two nodes at the ends of the segment needs to be measured to obtain the elastic energy. A schematic representation of a spring under compression is shown in figure 2.1.

While the elasticity of the fiber segments is important, the fibers also have a bending interaction. Resistance to bending, or bending rigidity, is common to biopolymer fibers and is one of the properties that can grant stability to a sub-marginal network. Contrarily to stretching, bending interactions depend on two consecutive segments, or three consecutive nodes, of the same fiber. The reason

Figure 2.1: A spring representing a fiber segment in rest (1) and under com-pression (2).

for this is that the bending energy is measured through the angle of deflection at the node connecting to fibers. This is shown in figure 2.2.

Figure 2.2: Two consecutive fiber segments. In (1), the deflection angle is zero and there is no bending energy. The bending energy in (2) is non-zero as there is a deflection angle δθ.

The Hamiltonian for linear bending interactions of one fiber at a single node is given by the bending energy from the wormlike chain (WLC) model [21], where the same assumptions about linearity are made as in the stretching energy. It has the following form:

Hbending = 1 2 κ l0(δθ) 2_{. (2.2)}

Here, κ denotes the bending rigidity, and δθ represents the angle the segment makes with respect to the rest angle as seen in figure 2.2. It is important to note that, even though multiple fibers can be connected at one node, bending energy is only calculated for consecutive segments of the same fiber, as calculating the bending energy between two segments of different fibers has no physical meaning. The last node, at the end of a fiber, has no bending energy.

The bending rigidity and the stretching modulus are material properties of the fiber. They depend on the fiber thickness, with thicker fibers being stiffer. The relation between κ and µ and the fiber thickness can be approximated by the following expression for rodlike segments:

κ= π 4r

4_{E (2.3)}

µ= πr2E, (2.4)

with r the radius and E the Young’s modulus of the rod [22] (see section 2.4.1 for more information on the Young’s modulus). In order to do realistic network simulations, accurate values need to be obtained from experiments.

We can combine µ and κ into a dimensionless ratio that can be used to characterize network behavior. This ratio, which is made dimensionless by mul-tiplying with a factor of l−2

0 , is defined as: ˜ κ= κ µ 1 l2 0 . (2.5)

It represents the relative importance of the bending length scale lb to the

geo-metric length scale l0. A value below 1 indicates that it is energetically favorable

for a fiber to bend in order to comply to deformation, while a value above 1 implies that stretching is preferable. This stems from the fact that fibers in the network relax to the lowest possible energy configuration, and thus the lowest energy interaction is preferred.

Because we set l0= 1 in the simulations, only the ratio of κ to µ is important.

To simplify this even further, in the simulations µ= 1 is used, so that κ is the only parameter in determining the ratio of the bending and stretching energy costs. Note that there are no thermal contributions to the Hamiltonian, because the fiber thickness is too big for thermal fluctuations to have any noticeable effect. Moreover, the strains exerted during compression and as a result of remodeling are much stronger than any caused by thermal contributions to the energy.

2.2

Macroscopic interactions: Connected fibers

Up until now, only single fiber interactions have been considered. While only two forms of interactions on the microscopic level are described above, the macro-scopic response of such networks can already exhibit rich behavior when de-formed. To get a better understanding of how simple rules at the microscopic level lead to nontrivial behavior at the network level, the arrangement of fibers into a network needs to be discussed. Since the goal of this study is to com-pare simulated networks to fibrin networks, the choices that have been made in constructing the network reflect this.

Before building a full network, the rules for connecting individual fibers must be studied. In biopolymers, there are different ways in which (parts of) fibers may be connected to each other. The first of these is the branching of the fiber bundle. This means that, at a branch point, one fiber is split into two, where there exists an angle between the two resulting fibers. Several branch points in fibrin networks are shown in the first panel of figure 2.3. The connectivity of a branch point is three.

In addition to branch points, two separate fibers can also form connections. The two fibers can be connected in different ways. For example, two fibers can cross-link with help of a cross-linking factor [8]. While the strength and type of connection may vary in real biopolymers, in this study there is only one type of cross-linking in building up the lattice; the cross-linking is simply represented by a node that connects segments of two different fibers. An example of cross-linking in a real biopolymer network, along with a schematic representation of cross-links and branch points, is shown in figure 2.3.

Since the connectivity at cross-links is 4 and the connectivity at branch points is 3, the average connectivity of biopolymer networks lies in the range of 3< z < 4. This number is difficult to accurately measure experimentally, and is approximately z≈ 3.4 for collagen [10].

Figure 2.3: (a) Scanning electron microscope image of a reconstituted fibrin network with highlighted branch points (red) and cross-links (yellow). The scalebar indicates 2 µm. EM image courtesy of Karin Jansen (FOM Institute AMOLF).(b) Schematic network representation with highlighted branch points and cross-links in the same color scheme. Note that not all fibers are connected to limit local connectivity to 4(see section 2.3.4).

2.3

Choosing the network structure

The first important choice in creating a network is how to define the layout of the network, regarding both geometry and topology. There are several alternatives for the layout: The layout can be based upon one of multiple lattice types, to which disorder is then later introduced; a network orientation can be generated randomly; or the complete layout is obtained through scans of a sample. While the latter choice is not as easily implemented as the other two, because it requires precise measurements and segmentation software, it can provide useful insights in the real structure of a fibrous network [23, 24].

Naturally, a network based on a lattice is much more regular than a purely randomly generated network, which is often represented by a mikado network [25, 8, 26]. In a mikado network, monodisperse fibers of arbitrary length are randomly brought together and placed inside a box. Hingable nodes are intro-duced at the places where fibers cross. It is unsurprising that it is easier to manipulate a lattice-based network, because it is much more regular: a single manipulation can be repeated multiple times throughout the lattice.

Another important parameter in determining what the network will look like is the number of dimensions; a 2D network will have a simpler structure than a 3D network. Additionally, adding a third dimension adds more degrees of freedom to the network, lengthening calculations. Fortunately, results obtained from 2D network simulations agree well with those of 3D networks, as shown by Licup et al. [27].

2.3.1

Maxwell’s criterion

Before elaborating on the choices made in creating the network, it is good to revisit Maxwell’s criterion. Maxwell stated that, for a network of Hookean springs that cannot bend and only hinge at their connecting nodes, the network is stable if the average connectivity satisfies z≥ 2⋅d, where d is the dimensionality

[9].

The origin of this criterion lies in the fact that in a network of Nd 0 nodes,

there exist d⋅ Nd

0 translational DOF. When the number of constraints in the

network is equal to the number of internal DOF, the network is marginally stable. Since each spring adds one constraint to the network, and each spring is connected to two nodes, the amount of connections needed to obtain marginal stability in the network is two times the number of nodes in the system: 2d⋅N₀d. Dividing the total amount of required connections by the total number of nodes yields the following result for the average connectivity:

z= 2d ⋅ N₀d/N₀d= 2d. (2.6)

As long as the network has an average connectivity equal to this thresh-old, the network is called isostatic; isostaticity implies that the network has the minimum number of constraints to provide resistance to deformation in any con-figuration. The resistance in this case comes from the fact that no configuration can be found that does not compress or extend any springs. If more constraints are added, the network is considered overconstrained. An overconstrained net-work is stable, and simply has an increased resistance to deformations with re-spect to a (sub-)isostatic network. Sub-isostatic networks are underconstrained and show floppy behavior, providing no resistance to deformations. In the sub-isostatic phase networks can, however, be stabilized by applying deformations. These deformations effectively reduce the number of degrees of freedom.

2.3.2

Creating the network from a lattice structure

The networks presented in this thesis are based on a 2D lattice structure. More precisely, the structure is quasi-three-dimensional, as will be explained in the following section. The lattice base that was chosen is the equilateral triangle. An advantage of using a triangular lattice is that it is relatively easy to bring the average connectivity to the desired range of 3 < z < 4, by allowing fibers to overlap at nodes without being connected; Another one is that the fiber segments in a triangular lattice line up to form longer fibers that are also found in biopolymer networks. Additionally, triangles have a finite resistance to small shear deformations, unlike squares, at least when applying small enough shear strains along one of the sides of the square. The following is illustrated in Figure 2.4.

The lattice is specified to have the same amount of nodes N0 in both the

x-direction (horizontal) and y-direction (vertical), causing a full lattice to have N = N2

0 nodes. Note that, due to the geometry of an equilateral triangle, the

width of the network does not have to be equal to its heigth. In this study, the lattices are oriented such that one side of the triangle is perpendicular with the x-axis. Therefore, the width of the lattice, LX, is equal to the number of nodes

in one direction, N0, times the segment length:

LX= N0⋅ l0. (2.7)

The heigth of the lattice, LY is given by the number of nodes multiplied by the

height of the equilateral triangle: LY = N0⋅

√ 3

Figure 2.4: Comparison between the resistance of a triangle and a square when applying identical (small) shear deformations. (a) Applying a shear deformation to a triangle requires work to stretch one side and compress the other (shown in red). There is no energy cost associated with this deformation for the third side, which is therefore shown in green. (b) The sides of a square are not strained for small shear deformations, and thus there is no resistance to small deformations (all sides are green). In this comparison, the small shear angle was exaggerated for improved clarity.

Of course, by this definition of LX and LY, the size of the lattice includes

segments that are seemingly unconnected, as N0consecutive nodes only require

N0−1 connecting segments. The explanation for this is that the segment that is

sticking out at the last node of each row is connected to another node through periodic boundary conditions. These will be discussed in section 2.3.5.

2.3.3

Midpoints

The nodes at the edges of the equilateral triangles that make up the network are not the only nodes present. To facilitate bending within triangles, there are nodes in the middle of the triangles’ sides. Therefore, these nodes are re-ferred to as midpoints in the rest of the text. Since midpoints do not count as a crossing, the average connectivity of the lattice z remains unaltered. The reason why these midpoints are important will be discussed in the section re-garding remodeling (section 2.6.2). Finally, when midpoints are introduced to the network, the fiber rest length l0should be halved, to accomodate to the new

distance between nodes. Meanwhile, the mesh size remains unaltered. Due to the increased number of DOF in a network with midpoints, such a network has a lower stiffness than a network without midpoints at small strains.

2.3.4

From a triangular lattice to a disordered network

Now that the choice has been made to focus on networks based on a triangular lattice, it is time to address the implementation of disorder to the network. Additionally, the average connectivity of a full triangular lattice, which is z= 6, needs to be brought down into a range of 3 < z < 4 to simulate the average connectivity in biopolymer networks. There are two methods to do this, and they are described below.

The foremost method of reducing the average connectivity is called dilution; in the process of dilution, fiber segments are randomly cut out from the network until the average connectivity of the network has been reached. Subsequently,

fiber segments that are only connected to other fibers at one end, the so called dangling ends, are removed. This is done because the dangling ends can always find a relaxed position within the network, and thus do not contribute to its physical properties, which would unnecessarily require computational power.

Figure 2.5: A schematic representation of a network based on a triangular lat-tice after dilution and phantomization. Phantomized connections can be iden-tified by the loop that is drawn around the node.

While randomly removing bonds until the desired value of z is ob-tained makes the network disordered, diluting too much may leave the net-work fragmented into unconnected clusters, especially for smaller net-works. This may affect the compar-ison with experimental observations. To prevent this, another method of reducing the average connectivity has been developed, one that does not re-quire the removal of fiber segments. This method is called phantomiza-tion, and adds a pseudo-three dimen-sional component to the lattice, by al-lowing fibers to cross at nodes, with-out being connected. This method randomly phantomizes one crossing fiber per node (out of three), discon-necting it from the other two at the

node. Phantomization instantly changes the average connectivity to z = 4, bringing it in the desired range of 3< z < 4, without the removal of any fibers. Finally, the network may still be diluted after phantomization. An example phantomized network is shown in figure 2.5.

2.3.5

Boundary conditions

Due to the finite size of the networks used in simulations, boundary effects can-not be ignored and need to be addressed. For this, periodic boundary conditions are introduced to the network. The boundary conditions in the network ensure that nodes that are at one boundary are connected to nodes at the other side’s boundary. This means that the nodes at the top boundary are connected to the nodes at the bottom and the left is connected to the right. Networks with these boundary conditions are thus effectively folded into a torus shape, as can be seen in figure 2.6.

Because the energy of each fiber segment is based on the position of the nodes, as explained in section 2.1, it is important that the distances between nodes are correctly measured at the boundaries. In an undeformed network, this is relatively straightforward, as the distance between two opposing boundary nodes is reduced by LX or LY; these are the lengths over which the network is

periodic in the x- and y-directions, respectively.

When applying deformations, above calculation of the distance between the periodically connected nodes does not apply. If, for example, the network is vertically compressed to a factor , the periodic distance LY should be modified

to LY,new= LY,old⋅ . The boundary conditions are also important in applying

Figure 2.6: Schematic representation of the boundary conditions. (a) Here, the distances between a node (full black circle) and three nodes (stars) at the periodic boundaries are indicated with an arrow. The corresponding three pro-jections (hollow circles) of the first node are shown for periodic connections along x (pink), y (red), and x+ y (green). (b) A schematic showing how the network is folded for these boundary conditions. (c) A 3D image of a torus (Image from:

https://upload.wikimedia.org/wikipedia/commons/thumb/c/c6/Simple Torus.svg/ 500px-Simple Torus.svg.png).

2.4

Extracting mechanical properties

With the network structure and the microscopic interactions defined, this sec-tion’s focus will be on the network’s mechanical properties. Before turning to the question how the mechanical properties are extracted in the simulated networks, it is useful to understand what information can be obtained from biopolymer networks during experiments. Therefore, this section will cover how the physical properties are extracted from both the real and simulated networks.

2.4.1

Measuring mechanical properties

The mechanical properties that are measured in experiments are used to define a network’s resistance to a particular form of deformation. In other words, if these mechanical properties are known, so is the behavior of the network in different circumstances. Before heading into how the mechanical properties are measured, it is useful to have a look at what exactly these properties are that define the behavior of the network.

In experiments, the network properties are measured by calculating the stress-strain relationships to the corresponding deformations. The mechanical property that is measured depends on the applied deformation and the bound-ary conditions. For example, by applying a small shear deformation γ one can measure the linear shear modulus, G. The linear shear modulus is obtained by dividing the shear stress (σs) required to apply this shear strain by the resulting

shear strain (γ):

G=σs

γ (2.9)

Note that this method is only valid for obtaining the linear shear modulus, as it divides the total shear stress by the total shear strain at any measurement point.

The values of other mechanical properties are measured similar to the shear case. If the perpendicular boundaries are relaxed, the resistance to a small

(uni-)axial compressing or extending deformation  is given by Young’s modulus: Y =σn

 (2.10)

Here, σn is the normal stress and  is the axial deformation applied. When

extended, most elastic materials shrink in the perpendicular direction(s) and vice-versa. The ratio of the resulting deformation in the perpendicular directions to the applied deformation is known as Poisson’s ratio [22]. In the case that the perpendicular boundaries are fixed or the Poisson’s ratio is zero, the Young’s modulus is also known as the linear elastic uniaxial modulus, U .

There are additional elastic properties that characterize the response of a material. However, it is possible to write them in terms of the others and they are not all independent. One distinct method of deformation that is still worth mentioning is bulk expansion and compression: when applying bulk deforma-tion, the same compression or extension is applied along all axes of the network. If the nonlinear response of a material needs to be measured, a different ap-proach is needed. This apap-proach relies on derivatives, and is used in determining the mechanical properties during simulations (see section 2.4.4).

2.4.2

The rheometer

The most common tool in establishing the bulk mechanical properties of biopoly-mer networks is the rheometer [14, 7, 4]. Rheometers allow measurement of stress-strain relationships for various different deformations. It functions through having two movable and rotatable plates attached to either side of a sample, while performing oscillatory motions. A schematic of this is shown in figure 2.7. As mentioned in the introduction, the samples are reconstituted from pu-rified fibers. However, the reconstituted fibrous networks are surrounded by water during measurements, effectively forming gels [7]. When conducting ex-periments, one has to account for the additional effects that come from the water that is present. This is important in establishing what boundary effects come into play - an experiment that allows water to freely flow into or out of the gel behaves differently than one that does not.

2.4.3

Energy and deformations in simulated networks

The main goal of simulating networks is to extract their mechanical properties, and compare them with those of (reconstituted) biopolymer networks. In order to do this, the differences and similarities in the analysis of mechanical properties for reconstituted and simulated networks need to be known. With the basis of the measurement of biopolymer networks already explained, the following will focus on obtaining the mechanical properties from simulated networks.

The deformations that are applied in this thesis’ simulations all rely on minimizing the mechanical energy of the networks, effectively finding a zero-temperature equilibrium position in the given state of deformation. This means, for example, that when stretching the network in one direction, the strain is spread out across all fiber segments that are aligned with the applied strain. Of course, the precise response differs for each network. However, the Hamiltonian for the network keeps the same form, which is the summation of the energies of

Figure 2.7: A schematic representation of a rheometer. A sample is placed between a bottom plate, and either a cone or another plate. The Peltier plate (bottom) is used to control the sample temperature, and is used to apply a rotational strain on the sample. The resulting torque is measured by the second plate. The distance between the two plates determines whether the sample is compressed or extended. During experiments, buffer fluid is allowed to flow into or out of the sample freely.

all microscopic interactions in the network (equations 2.1 and 2.2): Htotal = 1 2 ⎛ ⎝∑i,j µ l0(δl ij)2 + ∑ <ijk > κ l0(δθ ijk)2⎞ ⎠ (2.11)

Here, the summation indices i, j and k run over all nodes that are present in the network after dilution. However, only connected nodes, ij, contribute to the elastic energy, and only three consecutive nodes, ijk, on the same fiber contribute to the bending energy.

Whenever deformations are applied to these networks, the energy of the total network needs to be minimized to find the correct configuration. Energy minimization is an iterative process, during which the positions of the nodes in the network are constantly updated to obtain a lower total energy. Energy gradients that are calculated at each node are used in this process. This method is called conjugate gradient descent, the details of which fall outside the scope of this thesis [28]. However, it is good to note this gradient methods only finds the nearest local energy minimum. That is why we employ small intermediate deformation steps.

An insightful way of characterizing deformations for a two-dimensional net-work is a 2 by 2 deformation matrix. This matrix represents all transformations to the x- and y-axes of the network. An example matrix, for shear deformations along the x-axis, is shown below:

(xnew ynew) = ( 1 ±γ 0 1 ) ⋅ ( x y) = ( (1 ± γy) ⋅ x y ) (2.12)

Here, γ indicates the amount of shear deformation that is applied in terms of the shear strain. The displacement along the shear axis, x, is proportional to the heigth, y. The sign depends on the direction of the shear deformation.

The deformation described in equation 2.12 is an affine deformation. Affine deformations are uniform deformations throughout the network. Consequently, in a disordered network, the affine configuration may not coincide with the lowest network energy, whereas in full lattices it will, provided there is no fiber buckling. However, it is still useful to use affine deformations in disordered networks, as it can shorten the search for a network’s lowest energy configuration by providing an educated guess of the lowest energy configuration. Therefore, whenever a deformation is applied to the network, the current configuration of nodes undergoes an affine transformation before energy minimization.

It is possible to show different types of deformations in matrix form, even combinations of multiple separate deformations. However, it is important to realize that these matrices only depict the global (affine) deformations on these networks; the final node configurations inside the networks depend on their connectivity and microscopic interactions and are likely non-affine.

Actually, the application of deformations is done by acting on the network’s boundaries. Therefore, the periodic boundary conditions are an important tool in applying them. The periodic boundary conditions that are used in these networks are known as the Lees-Edwards (LE) boundary conditions [29]. LE boundary conditions apply deformation through adjusting the periodic distance between periodically connected nodes at the network’s boundaries, causing the whole network to adjust. More precisely, the periodic projections of the net-work are shifted based on the applied deformation. As an example, figure 2.8

Figure 2.8: A schematic representation of the LE boundary conditions in an (a) undeformed network and (b) sheared network. The projections (gray) in the sheared network have been shifted according to the shear deformation γ. This causes the distances along the periodic boundary to be consistent with the rest of the network. Example periodically connected node pairs are shown (full circle and stars). The distance between each pair is calculated with a projection (hollow circle). Note that the distances along the periodic boundary change with the deformation in (b).

shows how the LE periodic boundary conditions are used in applying shear deformations.

2.4.4

Mechanical properties in simulated networks

With the definitions of the energy and deformations, the mechanical properties can finally be extracted. Since the energy of the network is already calculated to find its configuration during deformation, a different approach for calculating the mechanical properties is used. The different properties, such as the linear shear modulus, G, can also be calculated by using energy derivatives. In fact, nonlinear stiffness can also be calculated with energy derivatives, which makes them useful when approaching the critical point when measuring nonlinear stiff-ness [16, 27].

The differential forms of the expressions for the moduli are related to their stress-strain definitions by the following expression for the energy:

E= ∫ ⃗F ⋅ d⃗s = A∫ σγ⋅ dγ. (2.13)

Combining above relation with the expression for the shear modulus yields the following result for the nonlinear shear modulus K:

K= dσs dγ = 1 A d2E dγ2. (2.14)

Note that for pure shear deformations the area, A, remains unaltered, while for other forms of deformations it may change, depending on the Poisson’s ratio.

Figure 2.9: Strain stiffening in different types of biopolymers and the synthetic polymer polyacrylamide. The biopolymers show nonlinear stiffening after a critical strain has been reached, whereas polyacrylamide does not show this effect. Image by Storm et al. [7].

2.5

Behavior of fibrous networks

Before describing the specific behavior of fibrin networks, it is good to have a look at more general fibrous network behavior. With help from the definitions for the mechanical properties, it becomes possible to understand the character-istic strain stiffening observed in fibrous networks.

2.5.1

Strain-stiffening in biopolymer networks

Biopolymer networks typically respond to strain by stiffening. This property prevents the tissues from being deformed, while still allowing for the flexibility needed. A tissue that is flexible for small levels of deformation is less susceptible to damage than a rigid tissue; structures that are too rigid are more prone to breaking under stress and strain [7]. A good example of strain-stiffening in the human body is arterial tissue expanding and contracting during a heart cycle, where the tissue is compliant to internal pressure at low vessel diameters and resistant at high vessel diameters [17].

A graph showing the stiffening of different biopolymer types under shear strain on a log-log scale is shown in figure 2.9 by Storm et al. [7]. Also visible in this figure is polyacrylamide, a synthetic polymer, which does not exhibit any stiffening.

The degree of strain-stiffening in athermal biopolymer networks is related to the values of the mechanical properties µ and κ. Usually, µ is much larger than κ [11], and the actual values of these depend on the fiber bundle thickness,

which depends on the number of constituent (proto-)fibrils. Expectedly, bending thicker fibers is more difficult than thinner fibers. However, thicker fibers are also more resistant to stretching and compression.

For small enough network deformations, the network can comply by bending fibers. As κ is orders of magnitude lower than µ, the energy cost, and thus the shear modulus, is orders of magnitude lower for smaller values of shear strain (see figure 2.9). While increasing the shear strain from zero, the network thus nonlinearly transitions from a bending-dominated towards a stretching-dominated stiffness regime. The cross-over phase is very interesting, as it shows a critical-like phase transition [16].

The critical-like phenomena that can be observed in these networks have been studied with experiments and simulations similar to those presented in this thesis, and show universal behavior for different network types; the data of these networks are overlapping after a Widom collapse [16], with parallels to the critical phase transitions observed in ferro-magnetic systems [30], where the temperature is replaced by the shear strain as a control parameter and the magnetization by the shear modulus. Even though the study of critical phase transitions has been a prominent research field within physics, it falls outside of the scope of this work.

2.5.2

Fibrin network stiffening

As mentioned previously, the motivation for this work lies in experiments carried out on reconstituted fibrin networks. During the experiments, fibrin exhibits stiffening after being subjected to (axial) compression. This stiffening is separate from the strain-stiffening in that it applies after the applied stress and strain are removed. Before showing the results from these measurements, the deformation steps that were applied to obtain these will be briefly covered.

The fibrin gels were reconstituted from isolated fibrin, and their mechani-cal properties were extracted by a rheometer. The basic principle behind the experiments was to slowly compress the sample up to a certain value of axial compression () in steps, after which the normal forces the sample exerted on the rheometer were allowed to relax. Performing these steps slowly is crucial in ensuring that the water that is inside the fibrin gel can freely stream in or out. This in turn affects the mechanical response of the gel: in case of slow compression, the water has ample time to flow out and only the mechanical properties of the remainder of the gel are measured. The resulting gel has a higher fibrin concentration and retains its size along the non-compressed axes (i.e. the surface area between the gel and rheometer plates does not change), behaving like a gel with a Poisson’s ratio of zero.

Measurements of the linear shear modulus were performed at each compres-sion step to measure the stiffening of the network in the gel. To see how large the contribution to the stiffening was for each value of , the sample was ex-tended back to = 0 after each compression step. At this stage, the linear shear modulus was measured again. The process of compressing and returning to zero was repeated multiple times, for increasing values of . Therefore, this process is called cyclic compression.

The results of the measurements of the fibrin stiffness during cyclic compres-sion are shown in figure 2.10. The stiffness is represented by the shear modulus, G. When initially subjected to compression, the fibrin gel shows a decreased

Figure 2.10: Lin-log plot of the shear modulus (G) in Pascals vs. the gap size (i.e. distance between rheometer plates) in µm of two reconstituted fibrin samples. Starting from a gap size of 1000 µm, the samples were compressed up to 80% during cyclic compression (black: multiple cycles; red: single cycle). Experiments by Bart Vos et al. (unpublished). The final stiffness for both samples is quite comparable, despite the number of compressive cycles.

shear modulus and becomes softer. This is due to the fact that the network is given more freedom to comply to the applied shear strain by bending, which is energetically favorable. After a compression step, the sample is brought back to zero extension, and an increase in stiffness that depends on the maximum amount of priorly applied compression can be seen.

What is especially interesting about this behavior is that fibrin does not show significant stiffening during compression, even if the gel has already been stiffened by a previous compressive cycle. This can be seen by comparing the shear modulus during compression for the single-cycle (red) and multi-cycle (black) data; for some reason, the stiffening can only be observed when the compression has been reverted.

A possible explanation for the stiffening behavior of fibrin is the following: fibrin networks make new connections during compression. When the new con-nections are subsequently stretched out as the gel is decompressed, they become stretched and increase the stiffness of the system.

In order to know whether forming new bonds (i.e. remodeling the network) can be responsible for the increased stiffness in fibrin networks, two questions need to be answered. First, it must be established whether fibrin fibers can form new bonds. Second, if fibrin can form new bonds: can the formation of new bonds explain network stiffening?

Figure 2.11: Schematic representation of an optical tweezer setup. The light of a highly focused laser beam that is pointed upwards is shown in red. The electric field is largest in the beam waist. However, the equilibrium posi-tion of the bead in the trap (trap center) is located slightly above the beam waist, because the photons falling on the dielectric bead transfer momentum onto it. The lateral displacement from the center of the beam gives rise to a force. For small displacements, this force is linear as for a spring, and so the force can be accurately measured through the displacement. Image from:

http://large.stanford.edu/courses/2007/ph210/sung2/images/f1.jpg.

2.6

Does remodeling explain fibrin network

stiff-ening?

Establishing whether fibers in a fibrin network can form new bonds has to be done via experiments. Since individual fibers are very small, with a diameter on the nanometer scale [1], the experimental method needs to be able to image individual fibers, as well as verify whether they can stick together. An exper-iment that meets these requirements was designed, based on optical tweezers. This experiment is described in the following section.

2.6.1

Probing fibrin fiber stickiness with optical tweezers

Optical tweezers are optical traps that can confine dielectric molecular beads in a highly focused laser beam with high precision. The dielectric beads are located in the waist (focus) of the laser beam, where the electric field is highest. Since the photons of the laser beam transfer momentum onto the bead, the center of the trap is actually slightly behind the beam waist. When displaced from the trap center, the beads experience a restoring force based on the magnitude of their displacement. The force exerted on the bead can be acurately measured through this displacement. An image that shows a schematical optical tweezers setup is shown in Figure 2.11.

Figure 2.12: Top-down fluorescence microscopy images obtained during the ex-periment are shown along with a schematic showing the current bead configura-tion. (a) Two pairs of two beads (labeled I + II and A + B) hold one fiber each at different heights. (b) The top fiber is rotated and (c) lowered to make con-tact. (d) The fiber attached to beads I and II is then moved horizontally, which also causes the other fiber to move. Finally, in (e) and (f) the effects on one fiber of moving the other in combinations of horizontal and vertical movements are shown. Experiment and images by Bart Vos et al. (unpublished).

For the experiment performed on fibrin, a four-trap optical tweezer setup was used in conjunction with fluorescence microscopy. The four beads were attached to the ends of two separate fibrin fibers that were obtained from a thrombin-fibrinogen solution. The fibers were then fluorescently labeled and attached to the beads.

To verify whether fibrin fibers can form new bonds between them, the two separate fibers were brought into a crossing position at different heights, by defocusing one pair of traps. The fibers were then slowly brought to the same height by focusing all four traps, until they were touching at the crossing. Sub-sequently, to see whether connections between the two fibers had been formed, one fiber was moved horizontally, expecting that if the fibers had indeed formed a bond, the other fiber would be forced to move as well. In addition to the fact that the other fiber also moved in the direction of the pull, simultane-ously performed force measurements confirmed that fibrin fibers can connect to each other. Fluorescence microscopy images obtained during the experiment are shown in figure 2.12, along with a schematic representation.

The strength of the newly formed bond has not yet been determined, since the binding strength of the bead to the fibers was too small to pull the two connected fibers apart: the fiber simply detached from the bead. However, it was already established that the fiber-fiber bond strength is at least in the order of hundreds of piconewtons.

2.6.2

Adding remodeling to the network simulations

The challenge in adding a remodeling mechanic to the simulations is getting a microscopic interaction to alter the macroscopic network behavior in the desired way. In addition to this, it is desirable to keep the total fiber volume constant during the remodeling process, to prevent an unphysical addition of material.

simple idea: if two nodes, that are part of different fibers, come within a certain distance d as a result of externally applied deformation, they should merge together and be counted as one node. In Figure 2.13, a schematic representation of this is shown. When using this manner of remodeling, the fiber volume is conserved, while the local connectivity is increased. Additionally, since nodes themselves do not represent any material, and the bending energy of different fibers can be calculated over the newly merged node, it is justified to remove one of the two merged nodes from the system.

Figure 2.13: Schematic of remodeling within one triangle. Two nodes approach each other to within merging distance due to compression (black arrows), which causes them to merge (red arrows).

Before the implementation of any kind of remodeling during which two nodes merge, the distribution of distances between nodes should be measured, in or-der to get a good estimate for the remodeling distance d. Since the distances between any combination of nodes should be measured, the function that calcu-lates distances should be expanded to also be able to calculate distances across the periodic boundary, and not only at the boundaries themselves. For the un-deformed network this is irrelevant, but especially at high values of compression this can give notable differences in the calculated distance.

To better explain how our extended LE periodic boundary conditions affects the distance calculations, a schematic is shown in figure 2.14. In this schematic, the difference in calculated distances for two nodes close to the boundary is shown. The old method calculates the distance in the wrong direction, causing it to overestimate the distance. Because the new method allows projecting non-boundary nodes, it can be ensured that the calculated distance is correct. The correct distance is simply the shortest distance between a node and either another node, or that node’s projection.

As mentioned before, the measured distance between nodes does not only affect the energy level and magnitude of the restoring force of the springlike segment connecting them; it also determines the direction. In extreme cases, miscalculating the direction can even lead to a network that (unphysically) crumples up. However, this problem only occurs in a network that allows for

Figure 2.14: Schematic showing how the distances between two nodes are cal-culated across the periodic boundary for the conventional and extended LE periodic boundary conditions. (a) When two nodes are at opposite boundaries, one node is projected with a shift of the size of the lattice (here: LY). In this

situation, both methods are identical. (b) When one or two of the nodes are not at opposite boundaries, the extended LE periodic (boundary) conditions can be used to calculate the correct distances by using an appropriate projection. (c) The original LE boundary conditions are not suitable for this, which can lead to an overestimation of the distance and obtaining the wrong distance.

Figure 2.15: The number of pairs within distance for different remodeling dis-tances d for a 20 by 20 network under compression. Higher values of d cause more node pairs to be within range of each other and a smaller compression onset before any pairs are counted.

remodeling, since nodes are connected across the periodic boundary only after the introduction of remodeling.

A useful way to characterize the remodeling distance d is expressing it in terms of the fiber rest length l0. In this way, a remodeling distance of d= 0.1

means that remodeling will take place if two fibers are within ten percent of l0

of each other. A simulation of a 20 x 20 network was performed in order to get the distribution of nodes within distance for d= 0.2, d = 0.1 and d = 0.05 during compression. The dimensionless bending-to-stretching-ratio was ˜κ= 10−3 and the average connectivity after dilution and phantomization was z≈ 3.4.

In figure 2.15, a plot of the number of nodes within different values of d are shown for an increasing value of compression . Of interest are the progression of the nodes within distance and the onset of the first nonzero value for each value of d. From experimental data, it is known that the stiffening of the fibrin networks occurs after a compression onset of  ≈ 0.1. Therefore, the value of d that will be implemented in the remodeling simulations should also have an onset of approximately ten percent compression. The value of d that best reflects this property is d= 0.05. Because the onset in larger networks might be slightly smaller than was obtained for a 20 x 20 network, a value in the range of 0.02< d < 0.05 was used in the simulations.

If the distribution of nodes within cross-linking distance d= 0.05 is split up into three interactions, based on whether the nodes involved are junctions or midpoints, e.g. midpoint-to-junction, it can be seen that different fibers are most likely to approach each other at midpoints (see Figure 2.16). It is thus expectable that if fibers could only form new connections between midpoints,

Figure 2.16: For distance d= 0.05, the number of pairs are plotted based on what type of node they are: midpoint-to-midpoint (blue), midpoint-to-nonmidpoint (black), nonmidpoint-to-nonmidpoint (red).

the network would still stiffen due to remodeling. An advantage of only allow-ing midpoint-to-midpoint mergallow-ing is that fibers are not allowed to connect at existing junctions, possibly leading to unrealisticly high local connectivities. To limit the maximum possible connectivity, no more than three nodes are allowed to merge into one.

Chapter 3

Results

In this chapter the results of the simulations are presented. The simulations have been performed on different 50 by 50 2D triangular lattices that have been phantomized and diluted to an average connectivity of z ≈ 3.4. The results are shown in a way that allows for easy comparison with experimental results of fibrin networks, as the model was developed to explain fibrin behavior and the experimental results are shown alongside. The choice for the dimensionless bending rigidity also reflects this, with our values of ˜κ= 10−3_{and ˜κ= 10}−4_often

used in simulations of biopolymer networks [27, 16].

For the simulations, the described experimental procedure of compressing the network before extending it back to its original size was followed. During these simulations, the shear modulus was measured after each applied defor-mation step to directly compare the stiffening with experimental results. The effect of the remodeling distance on the amount of network stiffening was also investigated; an increased d did not only increase the final stiffness, but also led to an increase in stiffening behavior during compression.

Since the remodeling gives rise to a normal force when the network is ex-tended back to = 0, it was investigated whether it was related to the shear mod-ulus. The simulation results indicate an apparent quadratic relation between σn and G. Experimental results are also indicative of a power-law dependence

like G− G0∝ σαn. The exponent was found to be α= 23 for both experiments

and simulations.

3.1

Network stiffening as a result of cyclic

de-formation

In figure 3.1, the average shear modulus of an ensemble of up to ten simu-lated fibrous networks with different randomization seeds during various levels of compression and decompression with remodeling distance d= 0.02 is shown for ˜κ= 10−3_{. Note that the average shear modulus does not have the same}

num-ber of contributions at every shown data point, since (unphysical) subzero shear modulus contributions were removed and not all simulations were completed for all ten seeds. The shown data was smoothed according to a moving average pro-cedure, during which data point were averaged with the previous and next two data points. There are several aspects of the stress-stiffening curves that are

worth mentioning. 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5

ǫ

G

ǫ

max=0.05

ǫ

max=0.10

ǫ

max=0.15

ǫ

max=0.20

ǫ

max=0.25

ǫ

max=0.30

ǫ

max=0.35

ǫ

max=0.40

ǫ

max=0.45

ǫ

max=0.50

Figure 3.1: Graph of the average shear modulus G as a function of the applied compression  for multiple networks, where the uncompressed networks were first compressed in steps of δ = 0.01 (crosses), before being returned to its original size in reverse steps (circles). Different colors indicate different values of maximum compression. The dimensionless bending rigidity was ˜κ= 10−3and the remodeling distance d= 0.02. Data points include an average of up to ten realizations of a network.

First and foremost, there is stiffening of the network after being compressed and decompressed. The final shear modulus is higher for networks that have been subjected to larger compressions, with an offset in stiffening of ≈ 0.1. The increase in stiffness during the decompressive phase appears to be nonlinear, with an increase in the stiffening rate towards the final steps of decompression. This can be attributed to the fact that remodeling can connect nodes that are in close proximity at high compressions, but far apart in an undeformed network. Therefore, when the remodeled networks approach = 0, the strains on the fiber segments of the remodeled nodes increase rapidly.

During the initial compression these networks all soften, since they networks have more freedom in complying with shear deformation by bending. Identi-cal findings have already been reported for unremodeled fibrous networks with similar values of ˜κ [31].

At the latter stages of compression with high max, the networks exhibit

stiff-ening behavior. This is due to their remodeling and has a similar argument for the nonlinear stiffening in the decompressive phase: the segments connecting the nodes that have remodeled at lower compressive deformations are increasingly strained when approaching high compressions.

For comparison, in figure 3.2 the shear modulus is plotted against the com-pression for an ensemble of networks with ˜κ= 10−4_{. Here, the data was not}

subjected to smoothing, as there was a higher number of data points available for the larger compression values (out of ten). Interestingly, both the initial softening and remodeling-induced stiffening in the compressive phase are much

less pronounced than for ˜κ = 10−3_{. While it is expectable to see less initial}

compressive softening for a network with ˜κ= 10−4_{, as it starts out softer than a}

network with ˜κ= 10−3_{, it is not obvious why it also shows less stiffening in the}

compressive phase. 0 0.1 0.2 0.3 0.4 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4

ǫ

G

ǫ

max=0.05

ǫ

max=0.10

ǫ

max=0.15

ǫ

max=0.20

ǫ

max=0.25

ǫ

max=0.30

ǫ

max=0.35

ǫ

max=0.40

ǫ

max=0.45

ǫ

max=0.50

Figure 3.2: Graph of the average shear modulus G as a function of the applied compression  for multiple networks, where the uncompressed networks were first compressed in steps of δ = 0.01 (crosses), before being returned to its original size in reverse steps (circles). Different colors indicate different values of maximum compression. The dimensionless bending rigidity was ˜κ= 10−4and the remodeling distance was d= 0.02. Data points include an average of up to ten realizations of a network.

3.1.1

The effect of the remodeling distance on network

stiffening

To ascertain that the stiffening of the network is caused by remodeling, and to demonstrate how the remodeling distance d plays a role in this, the stiffness during (de-)compression of one network is shown with four remodeling distances in figure 3.3.

As expected, the stiffening becomes more pronounced when the remodeling distance d is increased, since the probability of remodeling increases with it. However, it appears that the additional gain in stiffening is only present for the larger compressions (i.e. for higher max). This does not only reflect on the

final stiffness, but also on the stiffening at the end of the compressive phase, as can be seen in figure 3.3.

Figure 3.3: Shear modulus vs. compression during compression (crosses) and decompression (circles) for one network with different remodeling distances: (a) d = 0.02, (b) d = 0.04, (c) d = 0.06, and (d) d = 0.08. The remodeling-induced stiffening is larger for increasing values of d, both in the compressive and decompressive phases.

3.2

A power law relates network stiffness to the

normal stress

As seen in the previous sections, the remodeled networks are stiffer than before remodeling. By plotting the residual normal stress after compression and de-compression with respect to the normalized shear modulus increase for multiple values of maximum compression on a log-log scale, a power law dependence of the shear modulus on the normal stress was found (see figures 3.4 and 3.5).

The power law dependence can easily be indicated by a straight line, where the slope of the line represents the exponent b of the power law: G= G0+ χσnb,

where χ is an undetermined pre-factor. The lines in figures 3.4 and 3.5 suggest an exponent of b = 2₃ is possible. Despite the fact that the exponents seem consistent for these two values of ˜κ, there are insufficient data to rule out a slightly different exponent.

Since it is interesting to know whether this relationship is also present in experimental data, the experimental data on fibrin from Bart Vos et al. (un-published) were also analyzed for this. These data are shown in figure 3.6, alongside a line of slope −2₃, as in the previous figures.

−100 −10−1 −10−2 −10−3 100 101 102

Residual σ

n

(G

−

G

0

)/

G

0 Max. Compression 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

slope: −

23

Figure 3.4: Log-log plot of the normalized increase in shear modulus(G−G0)/G0

as a function of the residual normal stress σnafter different levels of compression

for an ensemble of network realizations with ˜κ= 10−3. A line with a slope of −2

3 is inset to indicate power-law dependence of G on σn. Only data points

with stiffness of G≥ 2G0are shown, so that only networks that have undergone

−100 −10−1 −10−2 −10−3 100 101 102

Residual σ

n

(G

−

G

0

)/

G

0 Max. Compression 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

slope: −

23

Figure 3.5: Log-log plot of the normalized increase in shear modulus(G−G0)/G0

vs. the residual normal stress σn after different levels of compression for an

ensemble of network realizations with ˜κ= 10−4. A line with slope −2₃ is inset to indicate a power-law dependence. The line also allows comparison with figure 3.4. Again, only data points with stiffness of G≥ 2G0 are shown, so that only

−10

0

−10

−1

−10

−2

10

0

10

1

F

n

(N)

(G

−

G

0

)/

G

0

slope: −

2₃

Figure 3.6: Normalized shear modulus increase (G − G0)/G0 vs. the residual

normal force Fn on a log-log scale for fibrin after compression and

decompres-sion. Through the data, a line of slope −2₃ is shown to indicate power law behavior with exponent b= 2₃. The experiments were performed by Bart Vos et al. (unpublished).

Chapter 4

Discussion and Conclusion

The results of the simulations performed with our model indicate that a straight-forward remodeling mechanism can be responsible for stiffening in fibrous net-works. Moreover, with proper choices for z and ˜κ, this model shows the same qualitative stiffening behavior as fibrin that is subjected to axial deformation, which makes remodeling the prime suspect for the stiffening in fibrin networks. Even though the model is capable of capturing the most important effect, it is not capable of capturing the full behavior of fibrin networks. Fibrin networks show stiffening based on the maximum compression they have been subjected to, regardless of the number of compression or decompression steps. This is in contrast to the simulated network, which can only stiffen due to the remodeling mechanic, since it is only capable of forming new bonds and not sever them.

As expected, we have seen an increased amount of stiffening for higher values of compression, as well as for the remodeling distance. While the strength of the remodeling in the simulations may be finetuned through d, it questionable whether it is useful to do a quantitative comparison between simulations and experiments, since the model was developed to merely explain the behavior of remodeling in fibrous networks. Nevertheless, it is remarkable that such a straightforward and general model is able to capture the behavior of a complex biological system.

An interesting result is the power law-relationship between the normal stress and the shear modulus after remodeling, with exponent b= 2₃. It is not yet clear what the underlying reason for this exponential relationship is, especially con-sidering that Licup et al.[27] have found a linear relationship between the normal stress and the shear modulus, although it has to be noted that they exclusively looked at shear deformations without remodeling. A possible explanation may lie in the fact that remodeling causes a normal stress along the x- as well as the y-boundary, and both normal stresses contribute. However, to investigate what the origin of this power-law is, it should be subjected to further study.