Mesoscale modelling of fibrin clots: the interplay between rheology and microstructure at the gel point

Abstract

This study presents a numerical model for incipient fibrin-clot formation that captures characteristic rheological and microstructural features of the clot at the gel point. Using a mesoscale-clustering framework, we evaluate the effect of gel concentration or gel volume fraction and branching on the fractal dimension, the gel time, and the viscoelastic properties of the clots. We show that variations in the gel concentration of our model can reproduce the effect of thrombin in the formation of fibrin clots. In particular, the model reproduces the fractal dimension's dependency on gel concentration and the trends in elasticity and gelation time with varying thrombin concentrations. This approach allows us to accurately recreate the gelation point of fibrin–thrombin gels, highlighting the intricate process of fibrin polymerization and gel network formation. This is critical for applications in the clinical and bioengineering fields where precise control over the gelation process is required.

---

1 Introduction

The viscoelastic properties of blood clots are biomarkers for blood coagulation in health and disease. They can reveal changes in clot structure, coagulation kinetics, and the processes of clot retraction and fibrinolysis.^1–4 The mechanical properties of blood clots are determined by the microstructure of the complex fibrin network, which is formed by the polymerization of fibrinogen into fibrin fibers. Clots with altered fibrin microstructure have varying susceptibility to fibrinolysis.^2,3,5,6 The fractal characteristics of incipient clots have been probed by oscillatory shear measurements of the gel point (GP) in whole blood,^2,7–9 revealing the characteristic interplay between the structural properties and mechanical response. The clot's mechanical response can change due to various conditions including venous thrombosis and coronary artery disease. In this context, the viscoelastic characterization of clots can be used in the screening and diagnosis of coagulopathies and in monitoring therapeutic interventions.

Fibrin–thrombin gels are a widely studied model system for blood clotting, providing insights into the clot's microstructure and mechanical properties.^3,8,10–12 Thrombin plays a crucial role in the coagulation cascade by converting fibrinogen to fibrin. Experimental studies have shown that thrombin concentration can significantly alter the rate and microstructure formation of fibrin gels.^8,10,13–15 Investigations using scanning electron microscopy (SEM) have shown that the concentration of thrombin significantly influences the characteristics of fibrin fibers, including their structure, length, and diameter. Specifically, lower thrombin concentrations result in thicker but less dense fibrin fibers, whereas higher concentrations lead to thinner, denser fibers. Similarly, the clots’ fractal dimension, D_f, has been identified to correlate with thrombin levels.^14,16–18

Chambon and Winter¹⁹ demonstrated that, at the gel point, both the storage modulus (G′) and the loss modulus (G′′) exhibit a power-law dependence on the angular frequency (ω), with the relationship G′′ ∼ G′ ∼ ω^α, where α represents a stress relaxation exponent. This relationship indicates that, at the GP, the phase angle (δ = tan⁻¹(G′′/G′)) becomes independent of frequency and equals δ = απ/2. Evans et al.²⁰ demonstrated that the GP of coagulating blood could provide a robust measure of a ‘clotting time’ and that power-law exponent (α) is sensitive to variations in thrombin availability. This sensitivity suggests that α can serve as a critical parameter for investigating and monitoring the impact of thrombin on the microstructure and stability of clots in a clinical setting.^2,20

Over the last decade, different numerical and experimental investigations have provided insights into blood clot rheological and microstructural characteristics.^12,14,18,21 In the context of fibrin clots at the gel point, relevant numerical studies have been also addressed. Curtis and coauthors^3,11 have explored the formation of incipient clots in fibrin–thrombin gels and heparinized blood using molecular dynamics simulations. They observed a correlation between gel-point time (t_GP) and fractal dimension, noting that D_f decreases as t_GP increases. Although their numerical model provides valuable insights into the microstructural properties of these gels, it does not account for variations of the viscoelastic properties with changing thrombin concentration. Overall, there is a noticeable knowledge gap in the numerical modeling of fibrin-clot formation. Existing models fail to simultaneously account for changes in the gel structure and viscoelastic response as thrombin concentration varies. Moreover, they do not adequately capture the key mechanical and structural relationships observed at the gel point. Here in, we provide a numerical scheme that captures relevant rheological and structural properties of fibrin–thrombin network formation. We demonstrate that the experimentally reported effects of thrombin concentration on the fractal dimension, gel elasticity, and gel time of fibrin–thrombin gels are accurately reproduced.^3,11

2 Definition of gelation system

Gel-network formation can be modeled using explicit coarse-grained representations of the aggregating molecules. These models are effective for capturing the specific interactions and kinetics of molecule aggregation,^22,23 as well as the distinct morphological features of these molecules.^24,25 However, the computational demands of such detailed simulations can limit their applicability to larger spatial and temporal scales. Here, we adopt a recently proposed mesoscale framework²⁶ to model biological clustering, accounting for hydrodynamic interactions and network formation kinetics. This approach balances detailed microstructural features and computational efficiency, enabling the study of morphology and rheological response of the gel over larger scales. The framework discretizes the gelling system into a particle-based representation, in which the microstructural elements are constituted by passive (P), active (A), and solvent (S) particles. The particle interactions are modeled using the smoothed dissipative particle dynamics method,²⁷ and the connectivity (bonding) between particles uses Morse potentials²⁸ (see simulation method for a detailed description of the model). As described by Zohravi et al.,²⁶ the rate of bond formation and maximum coordination number of bonds can be adjusted to model a variety of biological clustering mechanisms and morphologies.

Fibrin clot formation is a hierarchical and dynamic process. It begins with thrombin, generated through the coagulation cascade, converting soluble fibrinogen into insoluble fibrin. This leads to a lag phase where thrombin accumulates, followed by the polymerization of fibrin monomers into protofibrils and the cross-linking of fibrin fibers by activated factor XIII to form a stable, elastic network. The process involves dynamic interactions between fibrin monomers, which influence the clot's stiffness and contraction.^13,29 In Fig. 1(a), we provide a simplified sketch of the hierarchical fibrin-gel formation. It involves the aggregation of fibrin monomers into mesoscale-mesh structures, which subsequently polymerize into fibers and crosslink to form a network of fibers that constitute the clot.

Fig. 1 (a) Schematic of the fibrin-gel formation process, where fibrin monomers containing E domains (blue spheres) in the center of the monomers bind to D domains (red spheres) on adjacent monomers, forming initial mesoscale-mesh structures. Subsequent cross-linking between D domains stabilizes the fibrin mesh, creating a durable clot-network structure. (b) Particle representation for fibrin-network formation used in this study. Passive particles, P (dark blue), represent the initial mesoscale-mesh structures that can bond (blue bonds) to active particles, A (green), to mimic fibrin-network growth. Further bond (red bonds) formation between A particles models the cross-linking between fibrin fibers to consolidate the gel structure. Solvent particles (light blue) only interact hydrodynamically with the rest of the particles. (c) The initial concentration of passive particles (ϕ_int) serves as a proxy for thrombin concentration. At lower thrombin concentrations fewer mesoscale-mesh subunits form, leading to fewer fibrin fibers. Simulations are initialized with solvent and passive particles, and one active particle that acts as a seed to start the polymerization. A higher ϕ_int reflects an increased thrombin concentration, due to the larger probability of forming mesoscale-mesh subunits that can polymerize into fibers.

Here, P particles act as mesoscale-mesh subunits that polymerize into fibers by attaching to A particles. These active particles, once part of the fibers, can then cross-link to form the final gel structure (see Fig. 1(a)). Therefore, we can describe the fibrin-gel formation in two main stages: activation, where passive particles become active by bonding with active particles, and cross-linking, where these now-active particles bond among themselves to solidify the structure. Initially, only passive and solvent particles constitute the simulated system. The gelation initiates by placing a single active particle seed in the center of the domain.

Previous studies in fibrin–thrombin gels have shown that the fractal dimension, elasticity, and gelation time are sensitive to the activation and crosslinking mechanisms.^3,11,30 Based on the framework introduced by Zohravi et al.,²⁶ we investigate three representative clustering mechanisms, to identify the model that better captures the viscoelastic transitions during the gel formation dynamics. These mechanisms are denoted as branched, low-connected, and highly-connected, and allow us to sample different degrees of connectivity in the system and cross-linking rates.

During the activation stage, the bonding between passive and active particles occurs when particles are within a distance r* of each other. The bond formation process, repeats with nearby active particles, allowing for a maximum of three bonds per particle. During the cross-linking stage, the bonding between A particles also takes place when their distance is smaller than r*. Depending on the clustering mechanism a maximum number of “m” bonds can be formed per A particle. We set m = 0 for branched, m = 2 for low-connected, and m = 10 for highly-connected. This allows us to mimic the various degrees of connectivity between fibers.³¹

It is important to note that although activation and cross-linking can occur synchronously, they are consecutive stages from the P particle standpoint. As described by Zohravi et al.,²⁶ is it possible to control the rate of the different stages, using a lag time between activation and crosslinking, or equivalently changing the bonding probability of each stage. If the lag time is too short, the cross-linking stage can overtake the gelling process, leading to smaller aggregates that do not percolate the domain; conversely, if the lag time is too long, the formed gels get kinetically trapped in the branched condition. Here, we set a bonding probability of one for the activation stage – as long as the particles are within the bonding distance, they will form a bond –. Whereas, the bonding probability is reduced to 1/100 for cross-linking to mimic a slower cross-linking rate than activation. In the context of fibrin-gel formation, a lag time between activation and cross-linking can represent factors that lower thrombin levels or reduce the polymerization rate. These factors limit fibrinogen conversion and cross-linking.^32,33

We denote the initial concentration of P particles as ϕ_int = (N_p/N_t) × 100%, where N_p and N_t are the numbers of passive particles and total particles in the system, respectively. Here, we use ϕ_int as a proxy for thrombin concentration, and we investigate the effect of varying ϕ_int on the gelation process, as illustrated in Fig. 1(b). Thrombin is responsible for converting fibrinogen to fibrin, consequently generating the mesoscale-fibrin-mesh subunits. Therefore, it is expected that increasing thrombin concentration correlates with a higher number of initial passive particles. Overall, increasing ϕ_int implies that the probability of activation of passive particles increases. Similar approaches have been successfully used to model the effect of thrombin on the formation of fibrin–thrombin gels.^3,11 In the Results section, we show that ϕ_int effectively captures the structural and rheological effects of thrombin on the gelation process.

3 Simulation method

We use the smoothed dissipative particle dynamics (SDPD)^27,34 method to model the hydrodynamic interaction between the particles that constitute the system. SDPD is a thermodynamically consistent method that discretizes the fluctuating Navier–Stokes equations and has been widely used in the modeling of soft matter.^35–37 We consider a gelling system containing N particles with a volume V_i, such that , being d_i the number density of particles, r_ij = |r_i − r_j|, and W(r_ij,h) an interpolant kernel with finite support h and normalized to one. The evolution equations for the position of the particles is dr_i/dt = v_i, whereas the stochastic differential equation of the momentum, can be expressed aswhere v_i and p_i are velocity and pressure of the i-th particle. v_ij = v_i − v_j, e_ij = r_ij/|r_ij|, a and b are friction coefficients related to the shear η and bulk ζ viscosities of the fluid through a = (D + 2)η/D − ζ and b = (D + 2)(ζ + η/D). D is the dimension of the system. In eqn (1), we conveniently introduce the positive function F_ij = −∇W(r_ij,h)/r_ij. The last term in equation eqn (1), consistently incorporates thermal fluctuations in the momentum balance. SI.S1 and eqn (S1)–(S5) (ESI†) provide a detailed explanation of the terms A_ij and B_ij, the equation of state that defines the pressure p, and the form of the function F_ij.

The bonding between connected particles is modeled using a Morse potential²⁸ of the form U^bond_ij = D_Morse[1 − e^{−α(r_ij−r₀)}]² (2). The Morse potential is an ideal and typical anharmonic potential among the several molecular potentials, resulting in a bonding force F^bond_ij = −∂U_adh/∂r_ij = 2D_Morseα[e^{−2α(r_ij−r₀)} − e^{−α(r_ij−r₀)}] (3). Additional information about implementing the SDPD model, time window and dimensionless parameters used in the simulation are produced in SI.S3 (ESI†).

4 Gel characterization

The gelling mechanisms adopted allow us to model gels with different rheological responses and microstructural connectivity. To investigate the rheological response of the gels, we conduct the small-amplitude oscillatory shear (SAOS) analysis, as described in the following section. Whereas to account for variations in the microstructural features of the gels, we describe the connectivity in terms of their bond distribution, and the characteristic average number of bonds between particles- (N̄_B) and maximum-number (N_B-Max) of bonds per particle. In addition to the bond distribution, we also estimate the fractal dimension of the different gels, as an indicator of morphological distribution of the network. See SI.S2 (ESI†) for a detailed description of the fractal dimension computation. Given the statistical variability of the gels (and the stochastic nature of SDPD) depending on the initial distribution of the particles, the different gel characteristics are determined as mean values of five independent realizations of the gelling process.

4.1 SAOS analysis

We use SAOS^38,39 to study the linear viscoelastic response (complex modulus and loss tangent) for our model gels. We define the gel time (t_GP) as the time at which the sol–gel transition occurs due to network formation. The selection of the waveforms’ upper (5π 1/τ_SDPD) frequency limit is determined by considering the effects of fluid inertia. The lower limit, (0.25π 1/τ_SDPD) is selected for computational reasons. Lower frequencies pose numerical challenges, due to excessively long simulation times, or very small changes in the system behavior, making it difficult for our numerical method to detect meaningful responses. We first determine the ‘linear viscoelastic region’ for our gels by performing amplitude sweeps at a fixed frequency. In the linear viscoelastic limit, the moduli do not depend on the amplitude of oscillation. A strain amplitude of 0.05 was found to be well within the linear-viscoelastic response throughout the gelation process. To ensure a quasi-incompressible flow, the speed of sound was set to c = 50 ≫ V_max(ω) = γ₀ω for every frequency, such that the Mach number was always smaller than one. In addition, the solvent viscosity considered produces a Reynolds number Re = ρLV_max/η always smaller than one for all the frequencies investigated.⁴⁰ See SI.S4 and S5 (ESI†) for a detailed description of the simulation setup and boundary condition⁴¹ for stress measurement. To streamline the presentation of our results, we adopt non-dimensional time units, by normalizing the time with the corresponding SDPD viscous time, t_ν = (h)²ρ/η, where h is the kernel cutoff radius, ρ and η are the density and viscosity of the solvent, respectively. To account for the difference in bonding rate between activation and cross-linking, the lag time between the two stages corresponds to [small tau, Greek, macron] = 0.16. We chose to normalize the elastic modulus, viscous modulus, and frequency by utilizing the values Böhme and Stenger proposed in their work.⁴² As such Ḡ′ = ρ(s)²/2η²G′, Ḡ′′ = ρ(s)²/2η²G′′ and [small omega, Greek, macron] = ρ(s)²/2ηω are the dimensionless elastic modulus, viscous modulus and frequency. Whilst s is the gap between the fixed plate and the moving plate in the SAOS test. This facilitates the comparison with the results of different studies. By introducing the dimensionless frequency [small omega, Greek, macron] we account for variations in viscosity (η), density (ρ), and gap width (s). This approach allows [small omega, Greek, macron] to scale linearly with ω while adjusting for these changes in fluid properties. We conduct simulations for gel concentrations of ϕ_int = 7.5, 10, 12.5, 15%. Further details related to the maximum value of gel concentration can be found in Section SI.S2 (ESI†).

5 Results

We conduct an initial assessment of the connectivity and viscoelastic properties of the different gelation mechanisms (branched, low-connected, and highly-connected) to identify which model better captures the dynamics of fibrin-gel formation. In Fig. 2(a), we present the variation in the connectivity features for the three mechanisms at their final state. We define the final gel state as the point where bond formation stabilizes. Specifically, we consider the gel to be in its final state when the number of new bonds changes by less than 0.1% over 10 000 time steps. Their characteristic average bonds between particles and maximum number of bonds per particle can be summarized as follows:

Fig. 2 (a) Bond histogram of 3 proposed mechanisms. (b) Elastic (solid lines) and viscous (dashed lines) modulus of 3 proposed mechanisms at ϕ_int = 15%. Power-law slopes are included for the branched mechanism. The elastic modulus scales with a slope of two and the viscous modulus with a slope of one, suggesting a liquid-like behavior.

• Branched: N̄_B = 1.0 ± 0.01, N_B-Max = 3

• Low-connected: N̄_B = 2.0 ± 0.01, N_B-Max = 5

• Highly-connected: N̄_B = 5.4 ± 0.04, N_B-Max = 13

Regarding their rheological response, in Fig. 2(b), we compare the elastic and viscous modulus of the final gel obtained for each mechanism at ϕ_int = 15%. The branched mechanism shows a power-law rheological response with a low elastic modulus. The elastic modulus scales with a slope of two and the viscous modulus with a slope of one, suggesting a liquid-like behavior. Thus, despite its network-like connectivity this mechanism is insufficient to model viscoelastic response. In contrast, the plateau in G′ observed for low-connected and highly-connected mechanisms evidence a final higher rigidity of the gel. In particular, the highly-connected mechanism exhibits a larger change in the elasticity modulus during the gelation process. Besides the highlighted differences among mechanisms, in their bond count and elasticity modulus, the differences in their bond distribution reveal disparities in the rate at which the gel consolidates. For example, when comparing the highly-connected and low-connected mechanisms in (see Fig. 3) with the same average number of bonds, N̄_B = 1.86 ± 0.01, the highly-connected case exhibits faster gelation (t̄ = 0.63) compared to the low-connected (t̄ = 3.91) case. In this case, the larger number of bonds allowed during the cross-linking stage leads to a quicker consolidation or fast gelation of the highly-connected gel. Given that it depicts the solid phase in the highly-connected mechanism's final state, we can anticipate that it will be able to replicate the pre-gelation, gel-point, and post-gelation phases of a gelation process.

Fig. 3 Comparing the bond histogram for the highly-connected mechanism and the low-connected mechanism at the same average number of bonds N̄_B = 1.86 ± 0.01.

For further analysis and calibration of fibrin-gel formation, we focus on the highly-connected mechanism due to its ability to accurately represent the entire gelation process, including a wide range of elasticity modulus variations during gelling, and faster gel network consolidation. While the low-connected mechanism can model the gelation process, it is better suited for slower gelation processes with lower elasticity. It's important to highlight that the adopted gelling framework²⁶ supports modeling systems with features that span from low to high connectivity, by varying the degree of bond formation during the cross-linking stage. However, a detailed exploration of these parameters is beyond the scope of our current study. In subsequent sections, we address the identification of the three distinct gel phases—pre-gelation, gel-point, and post-gelation—using the highly-connected mechanism. Additionally, we compare the effect of gel concentration on the microstructural and rheological properties of the gel.

5.1 Gel phases identification

Identifying the gel point (GP) through viscoelastic characterization presents significant challenges. For accurate viscoelastic measurements, the degree of cross-linking in the gels should remain constant, or mutational effects should be negligible. The Mutation Number (N_μ) is a dimensionless quantity proposed by Winter et al.⁴³ It determines if a material's viscoelastic properties remain constant during the measurement. It specifically measures the rate at which material properties such as the storage modulus or loss modulus change in relation to the timescale of the rheological test. In our work, the Mutation Number is zero. However, gelling systems are inherently transient, changing as they gel. To address this, experimental methods such as halting the crosslinking reaction before conducting rheological measurements,^44,45 or employing multifrequency rheometric techniques, have been developed.^46,47 Numerically, a key advantage in identifying the GP is our ability to separate the gelling process from rheological characterization. Instead of having to pause or stop the reaction to monitor the GP during gelling, we first run simulations up to the post-gelation stage. We save the particle positions of the system at specific time intervals. Then, as a post-processing step, we perform SAOS tests to observe changes in the tan(δ) at various frequencies for selected samples from the entire gelling trajectory. We perform a SAOS test on one selected sample. If the system is in pre-gel (primarily liquid-like behavior), move forward in time to a later snapshot with higher N̄_B (indicating more bonds and progression toward gelation). If the system is post-gel (solid-like behavior), move backward to an earlier snapshot with a lower N̄_B. During the rheological characterization of these samples, we only consider existing hydrodynamic and bonding interactions, without allowing any new bond formations. This approach enables us to track changes in the system's viscoelastic properties at different stages, free from the ongoing gelling process. We essentially confirm that no structural changes have occurred by asserting that the Mutation Number is zero.⁴⁵

In Fig. 4(a) we show tan(δ) for the three phases of gelation – pre-gelation, gel-point, and post-gelation – at five different times (two at pre-gelation, two at post-gelation, and gel point) for the highest concentration, ϕ_int = 15%. In the pre-gelation phase, as the frequency (ω) increases, we observe a decrease in tan(δ), indicating behavior typical of a viscoelastic fluid. As gelation advances, the relationship between δ and ω weakens, and by the time we reach the gel point, tan(δ) becomes independent of frequency. Beyond the gel point, the material behaves as a viscoelastic solid, showing a characteristic dependency on frequency. In Fig. 4(b), Ḡ′ and Ḡ′′ evolution in these phases are shown. Both increase with time but G′ shows a more noticeable increment, indicating the development of the gel network. In Fig. 4(c), we show the values of D_f, N̄_B at the selected time sampled. Across the various systems evaluated, we identify that the average number of bonds between particles can be linked to the different phases of gelation. For instance, for ϕ_int = 15% a value of N̄_B < 1.86 ± 0.01 is indicative of the pre-gelation phase, while N̄_B > 1.86 ± 0.01 correlates with the post-gelation stages. The increase in G′ and G′′ and the decrease in D_f after GP indicates microstructural changes after GP. The GP occurs at lower values of both N̄_B and G′ and these parameters continue to change well beyond the GP. We must note that for a fixed ϕ_int, after the gel point the morphology of the network can keep changing due to cross-linking, without further increase of the gel volume fraction. This effect, leads to a morphological transition to more porous or open gels, along with a reduction in D_f.

Fig. 4 (a) tan(δ) at pre-gelation, gel-point, and post-gelation phases for ϕ_int = 15%, (b) Ḡ′ and Ḡ′′ evolution at pre-gelation, gel-point, and post-gelation phases for ϕ_int = 15% (c) D_f and N̄_B evaluation with time.

5.2 Effect of thrombin concentration on gel microstructure and viscoelasticity

We focus now on the variation of the microstructure and viscoelastic properties of the gel, with the initial concentration of P or equivalently to the thrombin concentration during the formation of fibrin–thrombin gels.

5.2.1 Microstructure behavior. Simulation snapshots of gel bond networks at GP are shown in Fig. 5(a) for ϕ_int = 7.5, 10, 12.5, 15% gel concentrations. The color map indicates the number of bonds per particle, evidencing a wider distribution in N_b for lower concentrations. These snapshots show higher concentrations lead to more densely packed structures and less open. Higher thrombin concentrations have been shown to cause the formation of dense microstructure and fibers with a high degree of branching due to rapid polymerization. These structures are often described using terms such as “fine,” “coarse,” “open,” “sparse,” and “tight”. Although our model does not explicitly model fiber branching at the monomer level, we can qualitatively say that higher initial concentrations of passive particles result in a denser network with more branch points due to the greater number of particles.^13,15 These effects on branching are consistent with SEM images of fibrin clots formed at different concentrations of thrombin, analyzed by Wolberg, as shown in Fig. 5(b).¹³

Fig. 5 (a) Snapshots of gel bond networks at the GP for ϕ_int = 7.5, 10, 12.5, 15% concentrations. (The color map indicates the number of bonds per particle for each gel concentration.) (b) SEM images of fibrin clots formed at different thrombin concentrations as imaged by Wolberg.‡¹³ (c) Mass scaling of the clusters as a function of their radius of gyration (SDPD units). We obtain the fractal dimension from these data D_f = 1.56 ± 0.1 for ϕ_int = 7.5% and D_f = 2.18 ± 0.04 for ϕ_int = 15% (d) D_f as a function of concentration compared with experimental results.¹¹ Dashed lines show the linear fitting of D_f, ≃0.2x + 1.8 for experimental and ≃1.2x + 1.1 for numerical results.

In general, a quantitative microstructure characterization is rarely performed.^8,23 In most cases, the results of quantitative characterization are limited to the average distance between branch points and the densities of the branch points. Here we provide a quantitative analysis of the microstructure complexity and branching of the gels by measuring their fractal dimension D_f. We compute D_f by considering the variation in the mass M of the formed clot with respect to their radius of gyration (R_g) as (see SI.S2.1, ESI† for a detailed description). The estimation of D_f is restricted on the length scales R_g < ε, where ε corresponds to the size of a blob in which ϕ_blob = ϕ_int. Beyond this scale, the network is considered to be homogeneous, and D_f = 3.

In Fig. 5(c) we present the measured D_f for the lowest and highest concentrations evaluated, indicating the range in D_f attained with our numerical model, where D_f = 1.56 ± 0.1 for ϕ_int = 7.5% and D_f = 2.18 ± 0.04 for ϕ_int = 15%. At the lowest concentrations, we have low-mass gels with the largest number of bonds observed at the gel point, N̄_B = 4.14 ± 0.02. This is consistent with percolation models, network theory, and polymer physics, where increasing bond connectivity leads to more compact networks with lower fractal dimension.^48,49 As the particle concentration increases, D_f also rises because the gel occupies more space. Consequently, fewer bonds per particle are required to reach a percolating cluster in a densely connected network.

Within the range of ϕ_int investigated, we observe a quasi-linear relationship of D_f with ϕ_int. This is consistent with experimental observations¹¹ on the dependency of fractal dimension at the gel point with thrombin concentrations, as depicted in Fig. 5(d). For comparison, in Fig. 5(d) we include the corresponding linear fittings of the measured and reported fractal dimension. In Fig. 5(d), we normalize both thrombin concentration (ϕ_thrombin) and ϕ_int to allow a proper comparison between the experimental and numerical data. The thrombin concentration in the experimental data is normalized against the concentration where a plateau in the gel-time is reached,¹¹ corresponding to ϕ* = 0.11 NIH per ml. For simulation results, we use ϕ_int|max = 15% as a representative limiting value. We must note that for larger particle concentrations (ϕ_int > 15%), the simulated gels exhibit D_f > 2. Therefore, we select the threshold for the limiting concentration to ensure consistency with experimental evidence for thrombin clots, where a fractal dimension much larger than 2 is unusual (see SI.S2, ESI† for further discussion about the selection of ϕ_int|max = 15%). This normalization provides a ϕ_int:ϕ_thrombin mapping, aligning the maximum gel concentration simulated, to the maximum effective thrombin concentration, beyond which the gelation kinetics is not affected.

It is crucial to point out that the SEM images in Fig. 5(b) show how fiber thickness changes with concentration. While our model does not explicitly incorporate fiber thickness, it effectively captures the network's key behavior by focusing on particle concentration, which governs branching and connectivity. In real fibers, thickness influences the fractal dimension: thicker fibers tend to form fewer and more widely spaced contact points, resulting in less dense, simpler networks with lower fractal dimensions. In contrast, thinner fibers, with their greater surface area for interactions, create denser, more interconnected networks with higher fractal dimensions. Our model implicitly reflects these effects by considering particle concentration, consistent with experimental observations showing that higher thrombin concentrations produce denser networks with more branching points and higher fractal dimensions, as evident in the SEM images. Also since the fractal dimension quantifies the spatial distribution of mass relative to the radius of gyration, branching has a more pronounced influence on changes in D_f compared to fiber thickness.

We must note that besides the fractal dimension, other relevant structural features in our model such as (i) the characteristic gel size and (ii) branch point density, can be used to correlate numerical and experimental evidence further. Whilst such mappings of other structural features are out of the scope of the present study, we identify a consistent trend among the simulation results. For instance, the formation of larger clusters with increasing the concentration of passive particles (ϕ_int), is consistent with the variations in the size of the fibrin gels observed experimentally (due to higher thrombin concentrations). A systematic mapping of the cluster size evolution between simulation and experiments can potentially offer a link between the characteristic size of the discrete particles in our model and the mesoscale-fibrin-mesh subunits. In experiments, the increase in thrombin concentration results in a denser network with more branching points. In our model, correlations of this branching density can be further explored by tuning the maximum number of bonds during cross-linking.

5.2.2 Viscoelastic behavior. In Fig. 6(a) and (b), we summarize the viscoelastic characterization of the gels with varying concentrations. Values of the elastic G′ and loss modulus G′′ at the gel point exhibit a power-law relationship with respect to frequency for all the concentrations evaluated, as depicted in Fig. 6(a). The corresponding values of α are also presented for comparison in Fig. 6(a). In Fig. 6(b), we show the variation of α = tan⁻¹(G′′/G′)/(π/2) for the different concentrations, as a function of frequency. We could observe that α is frequency independent and approximately equal to values of power-law exponent shown in Fig. 6(a), indicating that our gels are approximately at the gel point.¹⁹ We observe that small variations in ϕ_int (7.5% to 15%) cause significant changes in the microstructure, resulting in noticeable alterations in the viscoelastic modulus and gelation time. Both G′ and G′′ increase with concentration, but the greater change in G′ leads to a reduction in tan(δ).

Fig. 6 (a) Values of elastic modulus Ḡ′, and loss modulus Ḡ′′ as a function of frequency for ϕ_int = 7.5, 10, 12.5, 15% concentrations. (b) α = tan⁻¹(G′′/G′)/(π/2) as a function of frequency. (c) t_GP as a function of concentration compared with experimental results. Dashed lines indicate the curve-fitting t_GP ≃ 3.7 exp(−1.3x) for experimental and t_GP ≃ 3.3 exp(−1.2x) for numerical results. (d) Ḡ′ (measured at ω = 6.28 Hz) at the GP as a function of concentration compared with experimental results (all values normalized based on their respective values observed at the ϕ_int|max = 15% for present results and ϕ* = 0.11 NIH per ml for experimental results).¹¹ Dashed lines indicate the linear fitting G′ ≃ 1x + 0.1 for experimental and G′ ≃ 1.1x − 0.1 for numerical results.

For the simulated fibrin clots, we identify that the range of power-law exponent measured, α = 0.78–0.89, is in general larger than the ones reported for fibrin clots.^11,50 In our model lower power-law exponent can be obtained by (i) adjusting the fluid viscosity of SDPD particles, to reduce the viscous contributions – reducing G′′ –, or (ii) increasing the elasticity of the gel – G′ – with further fine tuning of the bonding potential adopted.

Overall, our simulations show that the changes in the viscoelastic properties of the gels, specifically the reduction in the power-law exponent and the rise in the storage modulus satisfactorily reproduce the response induced by thrombin in fibrin gels.¹¹ In Fig. 6(c) and (d), we further compare the values of t_GP and Ḡ′ (measured at ω = 6.28 Hz) as a function of ϕ_int, with previously reported evidence.¹¹ Following the normalization used rationalize the variation in D_f, in Fig. 6(c) and (d) the parameters t_GP and G′ are normalized using the values gel-point time and elasticity modulus observed at the thrombin concentration, ϕ* = 0.11 NIH per ml, corresponding to the limiting concentration where the gelation time reaches a plateau. For simulation results, we use ϕ_int|max = 15% as the limiting value. In Fig. 6(c), we identify that t_GP quickly decreases as the concentration of thrombin increases, indicating faster gelation at higher concentrations. The variation in gelation rates—from slow at ϕ_int = 7.5% to fast at ϕ_int = 15%—is attributed to differences in the activation probability of the highly-connected mechanism, which speeds up the cross-linking process. Similar to experimental results, the reduction in gelling time reaches a plateau, where further addition of thrombin does not shorten the gelling time. In this context, the limiting concentration, ϕ_int|max = 15% in our simulations is sufficiently high to capture the plateau in the gelation time (reflecting a stabilization of the gel formation), while maintaining the values of D_f < 2, as discussed in the previous section.

Comparing the functional dependence of t_GP, with the thrombin concentration, both experimental and numerical results show that at higher concentrations the t_GP can approximated by an exponential decay of the form t_GP ∼ ae^{−b[small phi, Greek, macron]}, where [small phi, Greek, macron] is the normalized concentration, a,b = 3.7,1.3 for experiments, and a,b = 3.3,1.2 for simulations. The corresponding fitted forms are included as dashed lines in Fig. 6(c). Experimentally, the existence of a limiting value of t_GP is attributed to the fact that in those experimental essays, thrombin concentration solely dictated the kinetics of fibrinogen activation. In this sense, the variation on ϕ_int in our model is equivalent to the variation in thrombin concentration in the experiments. At lower thrombin concentrations, the experimental evidence shows that the t_GP asymptotically increases. In our simulation results, we capture this trend for the range of ϕ_int evaluated. We must remark that for ϕ_int < 7.5%, we do not reach the gel point within the simulated computational time, consistent with such asymptotic behavior in t_GP.

In Fig. 6(d), we note that contrary to the gelation time, the storage modulus increases with thrombin concentration. This trend suggests that gels formed at higher thrombin concentrations have a greater ability to withstand elastic stress, remaining within the linear viscoelastic region. Experimentally, it has been observed that G′ plateaus, indicating that further increases in thrombin concentration do not affect the gel's storage modulus beyond a limiting concentration. This plateau is reached at the concentration of ϕ* = 0.11 NIH per ml, coinciding with the condition where gelation time also plateaus. In our numerical results, such plateau in G′ is absent, suggesting that the materials elastic properties can still evolve, even after the t_GP has stabilized. However, within the range of the mapping (ϕ_int:ϕ_thrombin) adopted, the simulated results are consistent with the experimental data showing a linear dependency of G′ with the thrombin concentration as shown by the dashed lines in Fig. 6(d). Over the concentration range mapped, the corresponding linear fitting leads to a slope of ≃1, for both experiments and simulations.

5.2.2.1 Relationship of power-law exponent α and fractal dimension. The relationship between the power-law exponent α and the fractal dimension D_f, provides insights into the interplay between the viscoelastic properties and the microstructure of gels. This relationship offers a method to indirectly estimate the structural properties of gels through rheological characterization.²

In the context of polymeric gels, the pioneering theoretical work of Muthukumar⁵¹ has been used to describe the correlation between the power-law exponent, α and the fractal dimension.^11,52 In contrast to percolation theory,^48,53 where unscreened hydrodynamic interactions are present, Muthukumar's work⁵¹ incorporated screening effects that can influence the network's morphology. In general, both screened and unscreened models serve as bounds, with physical systems potentially lying between them depending on the extent of screening.

In Fig. 7, we present our findings of D_f and α alongside the theoretical predictions (α = d/(D_f + 2), d = 3) for unscreened hydrodynamic interactions, discussed by Muthukumar.⁵¹ In Fig. 7, we can observe that the power-law exponent decreases as the fractal dimension increases. Overall, our simulation results align with the reported evidence for percolation theory, where α ranges from 1 to 0.7.⁵¹ Although the absolute values of α are slightly higher than those predicted by the theory, the observed slope is consistent with the behavior of gels with hydrodynamic interactions. We believe this effect is inherent to the SDPD methodology used to model the gels, as hydrodynamic interactions between the particles forming the gel are always present.

Fig. 7 α as a function D_f for theoretical results (unscreened theory),⁵¹ and present results.

We acknowledge that discrepancies may exist between fractal dimension measurements of the stress-bearing network and those of the overall density distribution. Specifically, density-based approaches (which reflects the structural network) can include additional features that do not necessarily contribute to the material's mechanical properties. Our reported fractal dimension corresponds to the structural network. Therefore, the discrepancies observed in our results may also arise from this approximation. Alternative strategies to determine the fractal dimension or estimate the fractal spectrum associated with the stress-bearing network are beyond the scope of this work but will be addressed in future publications. Despite the discrepancies observed between the theoretical predictions and our numerical results, the identified correlations between D_f and α indicate that the interdependence between the viscoelastic properties and the microstructure of fibrin gels can be effectively captured within the framework of our model.

Conclusion

This study introduces a mesoscopic model that simulates the formation of fibrin–thrombin gels, focusing on their rheological behavior and microstructural characteristics. Our model effectively replicates the complex fractal structure and viscoelastic properties observed in fibrin–thrombin gels at the gel point. One of the key aspects of our model is its ability to account for variations in the concentrations of thrombin into the structure of fibrin gels. The resulting gels have fibers with fewer branch points and larger pores at lower thrombin levels. Conversely, higher thrombin concentrations lead to denser fiber networks with more branch points and smaller pores. This variation in fiber structure directly influences the stiffness of the fibrin network, with higher thrombin levels generally leading to increased stiffness.

Furthermore, our model demonstrates that increasing the initial concentration of passive particles—akin to increasing thrombin levels—results in a decrease in the gelation time (t_GP) and the power-law exponent (α), while the fractal dimension (D_f) and elasticity modulus (G′) increase. These trends are consistent with both experimental and numerical data from previous studies, showing values of D_f ranging from 1.56 ± 0.1 to 2.18 ± 0.04. Notably, our results not only align qualitatively with existing evidence but also show good quantitative agreement. Additionally, our numerical model can capture the interplay between the viscoelastic properties and the microstructure of fibrin gels.

By accurately capturing the dynamics of fibrin polymerization and network formation, our model offers valuable insights for clinical and bioengineering applications requiring precise gelation control. This mesoscale clot model also sets the stage for future Lagrangian heterogeneous multiscale modelling⁵⁴ of clotting processes under physiological flow conditions.

Author contributions

N. M. and M. E. conceived and supervised the project. E. Z. designed the computational experiments and conducted the simulations and data processing. The original draft of the manuscript was written by E. Z. and N. M. and all authors contributed to the final version of the manuscript. N. M. developed the numerical implementation. All the authors discussed and analyzed the results. All authors approved the final version of the manuscript.

Data availability

The data that support the findings of this article have been included in the main manuscript and as part of the ESI† available in the online version of the paper. The data point of simulated gels are available in the fibrinGels repository, https://github.com/BCAM-CFD/fibrinGels.git. Experimental SEM images of fibrin gels in Fig. 5(b) from A. S. Wolberg, Blood Rev., 2007, 21, 131–142 were included with permission from the publisher.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research is supported by the Basque Business Development Agency through the “Mathematical Modeling Applied to Health” Project. Also by the Basque Government through the BERC 2022–2025 program and by the Ministry of Science and Innovation: BCAM Severo Ochoa accreditation CEX2021-001142-S/MICIN/AEI/10.13039/501100011033.

References

1. J. W. Weisel, J. Thromb. Haemostasis, 2007, 5, 116–124.
2. P. A. Evans, K. Hawkins, R. H. Morris, N. Thirumalai, R. Munro, L. Wakeman, M. J. Lawrence and P. R. Williams, Blood J. Am. Soc. Hematol., 2010, 116, 3341–3346.
3. D. Curtis, M. Brown, K. Hawkins, P. Evans, M. Lawrence, P. Rees and P. Williams, J. Non-Newtonian Fluid Mech., 2011, 166, 932–938.
4. S. N. Stanford, A. Sabra, L. D’Silva, M. Lawrence, R. H. Morris, S. Storton, M. R. Brown, V. Evans, K. Hawkins and P. R. Williams, et al. , BMC Neurol., 2015, 15, 1–8.
5. J. Collet, D. Park, C. Lesty, J. Soria, C. Soria, G. Montalescot and J. Weisel, Arterioscler., Thromb., Vasc. Biol., 2000, 20, 1354–1361.
6. M. Wygrecka, A. Birnhuber, B. Seeliger, L. Michalick, O. Pak, A.-S. Schultz, F. Schramm, M. Zacharias, G. Gorkiewicz and S. David, et al. , Blood Adv., 2022, 6, 1074–1087.
7. P. Evans, K. Hawkins, M. Lawrence and P. Williams, Biorheology, 2008, 130.
8. D. J. Curtis, P. R. Williams, N. Badiei, A. I. Campbell, K. Hawkins, P. A. Evans and M. R. Brown, Soft Matter, 2013, 9, 4883–4889.
9. N. Badiei, A. Sowedan, D. Curtis, M. Brown, M. Lawrence, A. Campbell, A. Sabra, P. Evans, J. Weisel and I. Chernysh, et al. , Clin. Hemorheol. Microcirc., 2015, 60, 451–464.
10. S. L. Rowe, S. Lee and J. P. Stegemann, Acta Biomater., 2007, 3, 59–67.
11. D. J. Curtis, PhD thesis, Swansea University, 2012.
12. R. I. Litvinov and J. W. Weisel, Matrix Biol., 2017, 60, 110–123.
13. A. S. Wolberg, Blood Rev., 2007, 21, 131–142.
14. A. Sabra, M. J. Lawrence, D. Curtis, K. Hawkins, P. R. Williams and P. A. Evans, Clin. Hemorheol. Microcirc., 2020, 74, 147–153.
15. R. A. Risman, H. A. Belcher, R. K. Ramanujam, J. W. Weisel, N. E. Hudson and V. Tutwiler, Biomolecules, 2024, 14, 230.
16. A. Takahashi, R. Kita, T. Shinozaki, K. Kubota and M. Kaibara, Colloid Polym. Sci., 2003, 281, 832–838.
17. R. Bateman, H. Leong, T. Podor, K. Hodgson, T. Kareco and K. Walley, Microsc. Microanal., 2005, 11, 1018–1019.
18. M. Brown, D. Curtis, P. Rees, H. Summers, K. Hawkins, P. Evans and P. Williams, Chaos, Solitons Fractals, 2012, 45, 1025–1032.
19. F. Chambon and H. H. Winter, J. Rheol., 1987, 31, 683–697.
20. P. A. Evans, M. Lawrence, R. K. Morris, N. Thirumalai, R. Munro, L. Wakeman, A. Beddel, P. R. Williams, M. Barrow and D. Curtis, et al. , Rheol. Acta, 2010, 49, 901–908.
21. A. Zakharov, M. Awan, A. Gopinath, S.-J. J. Lee, A. K. Ramasubramanian and K. Dasbiswas, Sci. Adv., 2024, 10, eadh1265.
22. N. Moreno, J. E. Perilla, C. M. Colina and M. Lísal, Mol. Phys., 2015, 113, 898–909.
23. S. Yesudasan, X. Wang and R. D. Averett, J. Mol. Model., 2018, 24, 1–14.
24. N. Moreno, S. P. Nunes and V. M. Calo, Comput. Phys. Commun., 2015, 196, 255–266.
25. S. Yesudasan, X. Wang and R. D. Averett, Biomech. Model. Mechanobiol., 2018, 17, 1389–1403.
26. E. Zohravi, N. Moreno and M. Ellero, Soft Matter, 2023, 19, 7399–7411.
27. P. Espanol and M. Revenga, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2003, 67, 026705.
28. P. M. Morse, Phys. Rev., 1929, 34, 57.
29. A. S. Wolberg and R. A. Campbell, Transfus. Apher. Sci., 2008, 38, 15–23.
30. T. C. Day, P. Márquez-Zacaras, P. Bravo, A. R. Pokhrel, K. A. MacGillivray, W. C. Ratcliff and P. J. Yunker, Biophys. Rev., 2022, 3, 021305.
31. M. H. Periayah, A. S. Halim and A. Z. M. Saad, Int. J. Hematol.-Oncol. Stem Cell Res., 2017, 11, 319.
32. S. Kattula, J. R. Byrnes and A. S. Wolberg, Arterioscler., Thromb., Vasc. Biol., 2017, 37, e13–e21.
33. N. E. Hudson, BioMed Res. Int., 2017, 2017, 2748340.
34. M. Ellero and P. Español, Appl. Math. Mech., 2018, 39, 103–124.
35. N. Moreno, P. Vignal, J. Li and V. M. Calo, Procedia Comput. Sci., 2013, 18, 2565–2574.
36. D. A. Fedosov, M. Peltomäki and G. Gompper, Soft Matter, 2014, 10, 4258–4267.
37. D. N. Simavilla and M. Ellero, J. Non-Newtonian Fluid Mech., 2022, 305, 104811.
38. R. B. Bird, R. C. Armstrong and O. Hassager, Dynamics of polymeric liquids, Vol. 1: Fluid mechanics, John Wiley and Sons Inc., New York, NY, 1987.
39. C. W. Macosko, Measurements and Applications, 1994.
40. A. Vázquez-Quesada, M. Ellero and P. Espanol, Microfluid. Nanofluid., 2012, 13, 249–260.
41. N. Moreno and M. Ellero, Phys. Fluids, 2021, 33, 012006.
42. G. Böhme and M. Stenger, J. Rheol., 1990, 34, 415–424.
43. H. H. Winter, P. Morganelli and F. Chambon, Macromolecules, 1988, 21, 532–535.
44. F. Chambon and H. H. Winter, Polym. Bull., 1985, 13, 499–503.
45. H. H. Winter and F. Chambon, J. Rheol., 1986, 30, 367–382.
46. M. Geri, B. Keshavarz, T. Divoux, C. Clasen, D. J. Curtis and G. H. McKinley, Phys. Rev. X, 2018, 8, 041042.
47. R. E. Hudson-Kershaw, M. Das, G. H. McKinley and D. J. Curtis, J. Non-Newtonian Fluid Mech., 2024, 105307.
48. P.-G. De Gennes, Scaling concepts in polymer physics, Cornell University Press, 1979.
49. D. Stauffer and A. Aharony, Introduction to percolation theory, Taylor & Francis, 2018.
50. K. Hawkins, P. A. Evans, M. Lawrence, D. Curtis, M. Davies and P. R. Williams, Rheol. Acta, 2010, 49, 891–900.
51. M. Muthukumar, Macromolecules, 1989, 22, 4656–4658.
52. C. Chidambareswarapattar, P. M. McCarver, H. Luo, H. Lu, C. Sotiriou-Leventis and N. Leventis, Chem. Mater., 2013, 25, 3205–3224.
53. A. Aharony and D. Stauffer, Introduction to percolation theory, Taylor & Francis, 2003.
54. N. Moreno and M. Ellero, J. Fluid Mech., 2023, 969, A2.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4sm01126k

‡ Reprinted from Blood Reviews 21, A. S. Wolberg, Thrombin Generation and Fibrin Clots, 131–142, Copyright 2007, with permission from Elsevier.

---