Fracture mechanics analysis of fibrin fibers using mesoscale and continuum level methods.

Computational models for simulating and predicting fibrin fiber fracture are important tools for studying bulk mechanical properties and mechanobiological response of fibrin networks in physiological conditions. In this work, we employed a new strategy to model the mechanical response of a single fibrin fiber using a collection of bundled protofibrils and modeled the time-dependent properties using discrete particle simulations. Using a systematic characterization of the parameters, this model can be used to mimic the elastic behavior of fibrin fibers accurately and also to simulate fibrin fiber fracture. In addition, a continuum model was modified and used to obtain the individual fibrin fiber fracture toughness properties. Using this model and the experimentally available fibrin mechanical properties, we predicted the range of fracture toughness (1 to k P a m ) values of a typical fibrin fiber of diameter 100 nm and its critical flaw size to rupture (~4 nm), both of which are not currently available in the literature. The models can be collectively used as a foundation for simulating the mechanical behavior of fibrin clots. Moreover, the tunable discrete mesoscopic model that was employed can be extended to simulate and estimate the mechanical properties of other biological or synthetic fibers.

Introduction

Blood clots are structural entities that are necessary for prevention of significant blood loss during injury. Clots also form during pathological states such as deep vein thrombosis, myocardial infarction, and strokes. In all cases, the clot must be mechanically robust in order to withstand both static and fluid dynamic mechanical forces. A network of fibrin fibers is responsible for providing the structural integrity necessary for a clot sustaining large forces and being able to resist fracture into fragments. The fibrin fibers are composed of fibrinogen, which is the main protein found in the blood at approximately 1.5–4 g/l for human blood plasma. Structurally, fibrinogen provides the rigidity for fibrin fibers, and ultimately blood clots, as it endures high forces and sustains low strains when a fibrin fiber is stretched in uniaxial tension [1]. Fibrinogen is a hexameric molecule that is composed of dual α, β, and γ molecules connected to a central E region with terminal D domains [2]. Fibrinogen is ~46 nm in length, and assembles into fibrin via a highly complex enzymatic reaction that initiates when thrombin cleaves fibrinopeptides A and/or B on the fibrinogen molecule. Subsequently, fibrin monomers polymerize linearly and assemble laterally into a half-staggered manner to form fibrin fibers [3,4]. These fibrin fibers are important structural entities for the stability of fibrin clots and must withstand forces under tension and shear [5]. Most importantly, the fibrin fibers themselves must not rupture under force in the vasculature. If fibrin fibers rupture, there is a potential for engendering a flaw within the fibrin clot matrix, thus causing a high risk of clot rupture. Clot ruptures are correlated with strokes, myocardial infarctions, and pulmonary embolisms and hence understanding this process from a mechanical standpoint is of great significance to the medical and scientific community [5-12]. In addition, it has been demonstrated that stretching of fibrin fibers may affect the lysis, or dissolution, behavior of blood clots [13]. This may be due to a mechanism where plasmin cleavage sites are inaccessible on fibrin fibers during stretching. Studies such as this necessitate further investigations that clarify how stretching, and ultimately fracture, affect the mechanical behavior of fibrin fibers.

In our previous work, we have developed some computational models [14,15] that can be used to predict the polymerization and formation of fibrin clots under various physiological conditions. However, these models need to be extended to the mesoscopic scale and modified and for individual fiber deformation f studies. At the fiber level, Hadi et al. [16] have developed a multiscale collagen model that necessitates understanding and computing fiber failure to predict macroscopic failure at the tissue level. Using results from this study and similar multiscale studies on collagen [17,18], one can generalize the framework for understanding fiber network behavior and failure in other extracellular matrix materials, such as fibrin. A micromechanics and mathematical model apposite for predicting the macroscopic mechanical behavior of keratinous fibers was developed by Puglisi et al. [19], where they verified the model results using experimental cyclic tensile tests. At the single fiber level, the biomedical engineering community would greatly benefit from computational studies that address the fracture mechanical behavior of fibrin fibers [6,20-23]. These studies are important not only from a mechanics standpoint, but also from a quantitative and energetics standpoint. The results from the current study will also elucidate the critical flaw size to fracture a fibrin fiber and physiological range of fracture toughness, and ultimately may shed light on the critical flaw size and energy necessary to fracture a fibrin clot. The current study addresses this issue by utilizing the fracture stress of fibrin fibers from several experimental studies [24-26] and employs this as an input for an existing finite element model for prediction of the critical flaw size necessary to rupture fibrin fibers. In addition, we present a discrete time-dependent method for modeling a fibrin fiber using bundled protofibrils, which was used to compute fracture and mechanical properties using various system parameters. The methodology utilized in this study is important for computational and experimental researchers studying fibrin fiber mechanics and disorders relating to thrombosis and hemostasis.

Discrete modeling of fibrin fibers

A mature fibrin fiber possesses a diameter on the order of 100 nm and has around 150 monomers per cross section [25]. The elastic modulus of fibrin fibers has been determined from Atomic Force Microscopy (AFM) as 3.9–15 MPa, depending on the type of cross linking of the fibers [27, 28]. The actual process of fibrin polymerization into fibers is a multiscale, multistep complex process, which necessitates the use of large-scale computer cluster resources. Consideration of fibrinogen conversion to fibrin monomers, its self-assembly, protofibril formation, lateral aggregation and finally formation of the fibrin fibers is a computationally expensive and sophisticated process [4]. To circumvent this, we omitted the nanoscale characteristics from the calculations and modeled the fibrin fibers as a composition of protofibril bundles, which allows for effective simulation of the mesoscale mechanics of fiber elongation and fracture. To achieve this, we developed a similar strategy to shape based coarse grain modeling as in Ref. [29]. The original procedure of the coarse graining procedure employs a self-organizing neural network to place the beads in the appropriate location. In this work, the beads were placed at the approximate center of mass in the cloud of atoms.

A typical protofibril consists of 15–20 monomers [30], which constitutes a length of 370 nm–460 nm. This protofibril can be modeled as a straight flexible chain with 20 beads connected to each other through a spring network. A fibrin fiber can be then modeled by stacking multiple protofibrils in an annular arrangement and also longitudinally. A schematic of such a fiber arrangement is shown in Fig. 1a, with the fiber being stretched on both sides by a force F. Also, in Fig. 1a, a full discrete fiber model with a length of 2 μm (red) constructed from the protofibril models is shown. An additional length of fiber, L = 450 nm, was prescribed at both ends of the fiber for applying forces and for reduction of the end local effects (blue). The protofibril model consisting of 20 beads is also shown in Fig. 1a. A magnified view of the cross section of the fiber model is shown in Fig. 1b and c, with a diameter of 100 nm. The fibrin fiber was composed of protofibrils by arranging them radially and circumferentially with a spacing of 7 nm. This spacing was estimated by considering the effective diameter of a fibrinogen molecule (10 nm) and increasing protofibril density towards the core of the fiber [25].

Fig. 1.

Mechanical modeling procedure for a single fibrin fiber (a) Schematic of the fibrin fiber (top) under uniaxial loading is shown with an applied force, F. The computational model was prescribed with length L + 2L. The end regions (blue) were used for pulling and the middle region (red) is the test fiber. The protofibrils were modeled with 20 beads (bottom). (b) Cross-sectional view of the fibrin fiber model, with diameter D. (c) Magnified view of the fiber surface shows inter-protofibril gaps. (e) The interaction of protofibrils to each other was modeled using (d) a Lennard Jones (LJ) potential and the interaction within the same protofibril was modeled using (f) a harmonic potential.

The distance between the beads was 23.68 nm and a gap of the same size was considered while the protofibrils were arranged longitudinally. While circumferentially arranging the protofibrils, a longitudinal pitch of 225 nm was considered for alternate protofibrils. This allowed simulation of the packing of the protofibrils to form an intertwined fibrin fiber (Fig. 1b). The elasticity of the protofibril was simulated using a harmonic potential between the beads as shown in Fig. 1f. The inter-protofibril interaction was simulated using a weak Lennard-Jones (LJ) potential (Fig. 1d). The * indicates the equation is in dimensionless form, consistent with the units section of the LAMMPS software manual [31]. The graphical representation of both potentials in action is shown in Fig. 1e. The LJ potential is responsible for maintaining the intact fiber and also simulates the friction between different layers in the event of fiber elongation or bending. A total number of 19,560 beads was modeled in the 2 μm fiber system.

The objective of this discrete protofibril model was to simulate the stretching of the fiber from both ends and eventually to estimate the break strain, elastic modulus, and fracture strength of the fibrin fiber. The bond potential and LJ potential parameters were prescribed in such a way to predict the results matching with the experimentally observed values. A major challenge that we faced with the model was matching the pulling rate of the fiber in the simulations to those observed in experiments. In experiments, fibers are pulled using an AFM tip at a rate of 320 nm/s [32]. This means, to achieve a break strain above 100% (typical values in experiments [27]) for a fiber with a length of 2 μm, 6.25 s are required. With a time step of integration of 30 fs, this would require 208 trillion computational steps. Our preliminary tests showed that with an optimum number of 6 CPUs, for a 19,560 bead fiber, simulations of 32.8 μs/day could be achieved. This would require a physical time of 522 years, and hence it was impractical to perform fiber stretching and dynamic simulations using pulling rates specified in experiments. Therefore, in this study we employed higher pulling rates than used in experiments and used those to estimate fracture and other mechanical properties. Because the pulling rates were much higher in the simulations, it also allowed for determination of the elastic mechanical properties, since it is well known that polymers behave more elastic at higher strain rates and more viscous at lower strain rates.

A protofibril was modeled using 20 beads and represents 20 fibrin monomers in this case. This corresponds to a mass value of 340,000 g/mol (molecular weight for human fibrinogen) for each bead in the simulation. The neighboring beads were set at a radius of 7 nm (radial and circumferential spacing of protofibrils) and the cutoff radius, r, for simulations was taken as 10 nm. The cutoff radius, r, was used to reduce computational cost. In essence, r was used to introduce a cutoff length in the LJ potential such that interactions between particles at this cutoff length were negligible. The length scale was non-dimensionalized by cutoff radius, rc = 10 nm−. Energy was considered using kT units, where k is the Boltzmann constant and T = 310 K is the temperature. The interaction between protofibrils was simulated using the LJ potential (Eq. (1)) and the beads were connected using a harmonic potential (Eq. (2)). The remaining variables, parameters, and other necessary values are provided in Table 1, in reduced units (non dimensionalized units or LJ units).

Table 1

List of simulation variables.

Parameter	Notation	Units	Value
LJ energy	εLJ	kBT	1 to 5.5
LJ zero potential distance	σLJ	rc	0.6236
Cutoff radius	rc	rc	10 nm
Spring constant	K	k_BT/σLJ2	6198.3
Bond distance	r0	rc	2.36
Time step	dt	τ=mσ^2/k_BT	0.0013
Stress	S	k_BT/σLJ3	–

The σ parameter was estimated by considering the equilibrium distance (the distance at which the potential reaches its minimum) as inter-protofibril spacing (0.7rc). Then . In Eqn. [2], r0 represents the equilibrium bond distance. The elastic modulus of the fibrin fiber was estimated as ranging from 3.9 to 15 MPa, based on prior experiments from the literature [6,26,33-35]. Based on this information, a stiffness of 0.0007 N/m (6198.3 in reduced units) corresponds to 4 MPa, for a fiber with 100 nm diameter and 2000 nm in length. This spring constant applies to the entire fiber and not for individual protofibrils. For protofibrils, this value was lower for individual springs. Since the exact value of the stiffness is not known at the protofibril scale, a scaling parameter (K) was multiplied by the bond potential (Eq. (3)). By tuning this parameter, this provided the overall flexibility of the fibrin fibers.

At this stage, the main unknown variables of the system were stretch (pulling) rate, ε , and (K). The computational models were created using MATLAB [36] and simulations were performed using LAMMPS [37]. The results were visualized using OVITO [38] and post processing was done using customized MATLAB scripts. In this first study, we examined the sensitivity of pulling rate on the mechanical properties of fibrin fibers.

Study 1: pulling rate (strain rate) sensitivity

Arbitrarily, the inter-protofibril interaction strength, ε, was set to 4kT and K was set to 0.25. The pulling rate of 1 m/s corresponds to 0.36 in LJ units. A set of 17 simulations was performed with pulling rates ranging from 0.05 to 0.5 (LJ units). The elastic modulus (E = σ/ε) and break strain (ε = 100(L − L0)/L0) of the fiber at various pulling rates were estimated and plotted in Fig. 2. Breaking (fracture) of the fibers was determined both visually and analytically when the inter protofibril distances becomes more than 2.36r. Physically, this means there is a gap at a particular cross section and the broken protofibrils are far away beyond r, such that any reattachment is highly impossible. The results show that the break strain varies with respect to pull rate; however, the elastic modulus shows a slight decline with increasing pull rate. The elastic modulus is in the range of 30–40 MPa, as opposed to the experimental findings of 3.9–15 MPa. This is due to the higher values of ε or K, which was investigated next. From the pull rate results, the break strain exhibits its peak value at a pulling rate of 0.36 (~1 m/s). The low break strains observed in these simulations in contrast to experiments are due to the higher pulling rates which are on the order of 106 times faster in simulations.

Fig. 2.

Effect of pulling rate on the break (fracture) strain and elastic modulus of the fibrin fibers.

Study 2: inter protofibril attraction effects

From the pulling rate sensitivity studies, we obtained a value of 0.36 for maximum straining of the fibrin fiber. Going forward, the pulling rate constant was set to 0.36 in the latter studies. To determine the influence of the inter-protofibril interaction strength on the fibrin mechanical properties, the K value was maintained at 0.25 and ε was varied from 1 to 5.5. With these values, simulations were conducted to estimate the break strain and elastic modulus of the fibrin fiber. The simulation results are valid in the range of ε from 1 to 5.5. Outside this interval, the simulations were unstable. Results of force vs. deflection from these studies are shown in Fig. 3a, where a stick-slip mechanical behavior was revealed. This behavior is due to the number of beads present in the protofibril as 20. Increasing the number of beads reduces the step behavior without affecting the slope of the curves. An additional study using 40 beads was conducted to verify this behavior. From the force deflection curves, the values up to 35 nm deflection were considered for estimating the slope, which allows for calculation of the stiffness and elastic modulus. The break strain was estimated when the fibers completely separated beyond a distance of 2.36r. Fig. 3a shows the force vs. deflection curve with increasing values of ε.

Fig. 3.

(a) Force vs. deflection curve of the fibrin fiber while the strength (ε) of the inter-protofibril interaction was varied. (b) Effect of ε on the elastic modulus and break strain.

From the results, it was determined that the break strain is not very sensitive to the protofibril attractive strength; however, the elastic modulus shows a proportional increase with protofibril attractive strength. Fig. 3b shows that the elastic modulus increases steadily from 2 MPa to 40 MPa when ε was changed from 1 to 5.5. At a value of ε = 1.5, the break strain exhibits its maximum value. In addition, the elastic modulus exhibits its maximum value at ε = 1.5 in the range observed in experiments.

Study 3: protofibril stiffness sensitivity

By maintaining ε = 1.5 and pull rate as 1 m/s (0.36), next we varied the K parameter of the protofibrils. 25 simulations with K = 1/scale, where the scale ranged from 1 to 15, were performed to understand the sensitivity of the mechanical properties on protofibril stiffness. The results from the simulations are shown in Fig. 4a and b. The results show that the increasing stiffness of individual protofibrils actually causes a decrease in the elastic modulus and break strain. When K was set to less than 0.5, the elastic modulus was in the range observed in experiments [6,26,33-35]. Also, the break strain demonstrated its higher values in this region. Until this phase of the study, out of pulling rate, K, and ε, only one parameter was varied at a time while maintaining a constant value for the remaining parameters. These studies have provided information of the extreme range of values that can be used to simulate the system to produce experimentally similar elastic modulus values. The next step was to understand the variation of the elastic modulus when two parameters were varied simultaneously.

Fig. 4.

Elastic modulus and break (fracture) strain results of the fibrin fiber by varying the scale factor of the protofibril stiffness. The same results in two different representations with (a) 1/K and (b) K are shown here.

Study 4: elastic modulus landscape

To explore the elastic modulus landscape, the pulling rate of the fibrin fibers was fixed at a constant value of 0.36. ε was then varied from 1 to 5.5 and K from 0.125 to 1.80. Fracture simulations were performed, and the results are shown in Fig. 5. The break strain (engineering strain) was estimated from Eq. (4), break (fracture) stress using Eq. (5), and elastic modulus from Eq. (6).

Fig. 5.

Results from 80 different configurations of the system with varying protofibril stiffness and LJ potential attractive strength. Variation of (a) fracture strain, (b) fracture stress and (c) elastic modulus plotted for stiffness scale changing from 0 to 1 and LJ potential strength varying from 1 to 5.5.

In Equations [4-6], L is the final length of fiber at fracture, L is the initial length of the fiber (2 μm) (2 μm), P P is the force acting on the fiber at the time of fracture, A is the initial cross-sectional area c of the fiber, and dσ/dε is the slope of the stress-strain curve.

The fracture strain results (Fig. 5a) show that a high stiffness scaling factor and low protofibril attractive strength d reduced the break strain values. Higher values of fracture strain w were observed when the stiffness scale was approximately 0.2 and ε = 4. The fracture stress is observed in Fig. 5b, with a maximum value of 350 kPa. These fracture stress values are much lower than the experimentally observed values of 7–11 MPa [39], likely due to an oversimplification of the model from Equation [5]. Also, the per monomer peak break force is estimated as 8.37 pN per monomer at ε = 5 and K = 0.5. The elastic modulus is shown in Fig. 5c, and displays a defined trend with ε and K. For ε values less than 2, the elastic modulus was poorly estimated. The results show that the experimental values of 3.9 MPa–15 MPa reside in the light blue shaded regions of Fig. 5c (refer to colorbar for exact representation).

The results from these four studies indicate that with a pulling rate of 1 m/s, ε= 3, and a protofibril stiffness of 0.00035 N/m, we can simulate the physiological values of the elastic modulus of a fibrin fiber.

Continuum level estimation of fibrin fiber fracture

The fibrin fiber at the micrometer scale can be considered as a flexible cylindrical rod. An approximate estimation of the fracture properties of fibrin fibers was computed based on the 3D finite element solution of pipes and rods [40], which provides reasonably accurate stress intensity factors for cylindrical rods and pipes with circumferential elliptical cracks. Using these validated theories as the basis, we examined the fracture properties of a fibrin fiber subjected to uniaxial tensile loading that exhibit mode I deformation.

Consider a cylindrical rod subjected to an external stress S, with a circumferential crack as shown in Fig. 6. The mode I stress intensity factor (K) at the edge of the crack is given as:

Fig. 6.

Schematic of Mode I tensile loading in a fibrin fiber with an elliptical crack (left) and magnified cross sectional top view of crack geometry (right).

Here, S is the remote uniform tension stress, a is the depth of surface crack, Q is the shape factor of the crack, and F is the stress intensity boundary correction factor.

The shape factor Q is given by:

Here, c is half of the circumferential growth of the crack. The results show that the stress peaks when a/c = 0.6 and a/D0 = 0.35, given Q = 1.63 according to Raju et al. [40]. Here, D is the diameter of the fiber.

The stress intensity boundary correction factor (F) was estimated as a function of a/c and a/D. For a/c = 0.6, FF is given as:

Let us consider S, the fracture stress of the fiber with a diameter D. The length scale was non-dimensionalized as a* = a/D0 and the stress was normalized as S* = S/S0. Equation (1) becomes non-dimensionalized as , where . When S* > 1, the fiber fails under fracture due to the applied stress exceeding the fracture stress. The variation of failure stress for varying fracture toughness and crack depth was estimated using MATLAB [36] and shown in Fig. 7a, where the failure region is highlighted. Similarly, the fracture toughness was estimated for variations in S* and a* (Fig. 7b). The line shown in Fig. 7b shows the critical stress location, above which the fiber fails under fracture.

Fig. 7.

(a) The relationship between the crack size and fracture toughness. The contours show the ratio of induced stress to fracture (failure) stress.

From experiments, the elastic modulus of the fibers is found ranging from 3.9 to 15 MPa depending on the type of cross-linking of the fibers [27,28]. Also, the fracture stress S of the fibers was estimated from experiments, varying from 7 to 11 MPa [39]. Based on the literature, we considered a fibrin fiber with an average diameter (D) of 100 nm [25]. The variation of fracture toughness (K) was estimated by varying the crack depth (a) from 0 to 100 nm and fracture stress from 7 to 11 MPa, and the results are shown in Fig. 8. As expected, the results show that increasing fracture stress leads to increasing fracture toughness. The critical flaw size (minimum crack depth for fracture of a fibrin fiber) at a fracture toughness of was 4 nm (~4 fibrin monomers). Our results also suggest that for a typical fibrin fiber, the fracture toughness ranges from to .

Fig. 8.

Estimation of the fracture toughness of a fibrin fiber with varying crack depth and fracture stress.

The analytical solutions from these studies can be useful in developing continuum-based fracture models, which can be used for discrete simulations of fibrous clots and can be used to estimate its collective mechanical properties and understand the behavior under stress.

Discussion on current limitations and future improvements

The methodology used in this research study to develop the discrete model closely resembles the shape based coarse graining technique. However, instead of utilizing self-aligning neural networks, we employed the center of mass to place the beads with an equivalent mass. This model, in combination with the Lennard-Jones potential for inter-bead interaction and Harmonic potential for intra-bead interaction was validated by matching the elastic fibrin properties from Atomic Force Microscopy (AFM) experiments [24]. These interaction parameters were further optimized to correlate with the experimental results of the fibrin fiber properties.

In the present work, the mesoscale fibers were modeled as stratified composites, which possess a different complex arrangement. This can be improved by conducting studies with random arrangements and also staggering them at different pitch lengths. In the present study we used the LJ potential. A detailed study can be performed in the future with various types of potentials such as the Morse potential, many body dissipative particle dynamics force field, and a combination of Coulomb and other potentials. The mesoscale fiber model can be used for studies where more insights about microscale fracture of fibrin fibers and fibrin networks are to be studied. The tunable nature of this model will be highly beneficial for incorporating force field mapping from molecular dynamics and coarse grain simulations. This can improve accuracy of the mesoscale model and can be employed to predict individual fibrin fiber fracture as a constituent of fibrin networks.

Our continuum model considers linear elastic fracture mechanics, which can be improved by accounting for viscoelastic properties and pulling rate sensitive parameters. The current continuum fracture modeling of the fibrin fiber was conducted assuming a linear elastic material, while in reality fibrin fibers exhibit nonlinear viscoelastic behavior. The model is assumed to be valid for the strain rates considered in this study, which were on the order of 106 times faster than those used in experimental AFM studies. This assumption can be improved by developing a model for fibrin fiber fracture based on viscoelastic energy formulations. The primary objective of the continuum model was to model the fibrin fiber at large continuum micro to macro length scales to correlate the structural fiber properties with experiments. This model is computationally efficient when compared with the mesoscale model, although there are no atomic and coarse grain details.

The continuum model was independently developed based on empirical formulation and is not correlated directly with the mesoscale model. Their connection can be modeled based on force matching or strain matching, which will be an improvement in predicting fibrin fiber fracture. In addition, a damage mechanics model could be used to study the fracture mechanical behavior of fibrin fibers, such as the one in Ref. [41], as a potential optimum solution.

Conclusion

In this paper, we have presented two different strategies for modeling the uniaxial tensile behavior of fibrin fibers and examining fracture mechanical properties. The first was a mesoscale protofibril discrete particle method and the second employed the utility of continuum methods. Using the discrete protofibril model, we systematically characterized the interaction parameters and accurately simulated the experimentally observed elastic modulus of a fibrin fiber. The model was also used to capture the rate-dependent fracture behavior of the fibrin fibers. Using the continuum model, we estimated the possible range of fracture toughness values of a fibrin fiber based on experimentally observed fracture stress and elastic modulus values. Our studies were used to determine the expected range of fracture toughness values of a typical fibrin fiber with an average diameter of 100 nm. The range was found to be between (1 to ) and the critical flaw size to rupture ~4 nm. To our knowledge, this is the first study to combine mesoscale and continuum level methods to determine the mechanical properties of fibrin fibers. These models can be modified to simulate the mechanical behavior of other biological fibers (i.e. collagen, elastin and fibronectin) and can be utilized to study the time-dependent fracture mechanism using computer simulations.