Upon exposure to biochemical stimuli that are either soluble species or cell surface-linked, platelet signal transduction results in calcium mobilization, which is a marker of platelet activation. Our model accounts for platelet activation via a neural network (NN) model trained in previous work using multicomponent agonist exposure data to determine a unique patient-specific, intra-platelet calcium mobilization response [19, 20]. The NN predicts the intracellular calcium concentration for any platelet, C a 2 + i t , based on the feedback vector (containing information from 1, 2, 4, 8, 16, 32, 64, and 128 s prior to the current instant) and the current concentration input of six agonists, namely ADP, TXA₂, collagen, thrombin, nitric oxide donor GSNO, and prostacyclin analog iloprost. The intracellular calcium concentration is then employed to determine the extent of integrin activation and adhesiveness of each platelet, ξ i t , given by ξ i t = ∫ 0 t C a 2 + i t − C a 2 + 0 d t , (1) where C a 2 + 0 is the basal level intracellular calcium ion concentration (100 nM) within a platelet. Detailed information on the NN architecture, training, and validation is available in Supplementary Section S1.

The Hill function is used to normalize the cumulative and recent-history calcium integrals between the basal and maximal levels of activation, α _min and α _max , to determine a time-dependent extent of inside-out signaling and integrin activation, F, given by F ξ i = α min + α max − α min ξ i n ξ i n + ξ 50 n , (2) where n is the Hill coefficient and ξ 50 is the critical value of ξ i for 50% activation.

Since NN computations of intracellular calcium concentration need to be carried out for each individual platelet within the simulation, we employ parallel particle decomposition to efficiently parallelize the NN module (see Section 2.2 for more details on the particle decomposition algorithm).

To track individual platelets as they move across the domain, aggregate to, or detach from the blood clot, we use the lattice kinetic Monte Carlo (LKMC) method [7–9]. The input for a LKMC simulation is a rate database for all possible events in the system at any given time. In the present case, the possible events include platelet motion (in the ± x, ± y, or ± z directions), bonding, and detachment. The rate of platelet motion in direction i is given by Ref. [21] Γ m o t i o n , i = D p l a t e l e t h L K M C 2 + max 0 , v ⋅ e i h L K M C , (3) where D p l a t e l e t , h L K M C , and v are the platelet diffusion coefficient in blood, the LKMC lattice spacing (1.5 μm), and local blood flow velocity vector, respectively.

The rates of platelet bonding/detachment are directly/inversely proportional to the extent of inside-out signaling, F (defined in Eq. 2), respectively. Furthermore, the role of the local shear rate in the shear-dependent breakage of ligand-receptor bonds and the effect of von Willebrand factor (VWF) are accounted for in the bonding and detachment rates. Detailed expressions for the bonding and detachment rates are available in Ref. [9].

With the specified rate database at system time t, the time step of the next event is chosen as Δ t L K M C = − ln ⁡ u Γ t o t a l , (4) where Γ t o t a l is the total rate of all events currently accessible in the system and u is a random number drawn from the uniform distribution in the interval (0, 1). The probability P i that event i with rate Γ i will be the next event is given by P i = Γ i Γ t o t a l . (5)

While there exist parallel open-source LKMC implementations in the literature, their applicability to model particle motion and aggregation in complex geometries is limited either by restrictions placed on domain geometries, or by limitations on the types of events that can be simulated [22, 23]. More importantly, they employ spatial domain decomposition, which would not be an efficient parallelization strategy for the present simulations. The reason for this is that the distribution of events, i.e., platelet motion, bonding, and detachment, with high propensities is extremely spatially heterogeneous; motion events are directly proportional to the local velocity field which is highly variable across the domain (see Eq. 3), and bonding/detachment events are restricted to a small portion of the domain. Since spatial decomposition parallel algorithms for LKMC involve cores performing events in their respective subdomains, and the LKMC time step is inversely proportional to the aggregate rate, a single core can end up being a bottleneck because time is advanced much slowly on that core [24]. Therefore, to efficiently parallelize the LKMC method with spatially heterogeneous (and time varying) event distributions, we have developed a custom in-house C++ code based on the message passing interface (MPI) library for inter-core communication that employs parallel particle decomposition instead of spatial domain decomposition. Here, each core is assigned a subset of particles, irrespective of their spatial locations, such that the sum of event rates is approximately uniform across cores. Moreover, to maintain this balance over time, new particles entering the simulation domain are assigned randomly among the cores. A schematic that illustrates the advantages of particle-based decomposition over domain decomposition is shown in Figure 2.

We have adapted the synchronous parallel kinetic Monte Carlo method developed by Martinez et al. for our parallel LKMC module [25, 26]. In the parallel algorithm, we first divide the system by distributing the particles equally among K cores. Then, the aggregate rate of all possible events on each core k , Γ t o t a l , k is calculated, and the highest aggregate rate Γ max is communicated to all cores, where Γ max ≥ max k = 1 , … , K Γ t o t a l , k . (6) A null event is assigned to each core with a rate Γ 0 , k given by the difference between the highest aggregate rate and the total rate on that core, i.e., Γ 0 , k = Γ max − Γ t o t a l , k . (7)

A null event is a “do-nothing” event where the system configuration remains unchanged. The introduction of these null events helps achieve time synchronicity, since each core now has the same total rate ( Γ max ) that determines the time step associated with an LKMC event. Each core now executes an LKMC event in parallel and independently of other cores, with time step and event probabilities according to Eqs 4, 5, respectively.

It can be readily seen that the parallel efficiency in the case of domain decomposition would be suboptimal due to the high fraction of null events in most cores (see also Figure 2), while parallel particle decomposition is designed to provide excellent load balancing. While load balancing is readily enforceable in particle-based decomposition, the chief drawback of this approach is that communication across cores becomes more challenging to execute because of the lack of well-defined inter-core boundaries. Ideally, any two particles assigned to different cores must not occupy the same lattice site, or perform lattice hops to or bonding events at the same lattice site. Strictly ensuring this condition after every event would require global communication across all cores after each LKMC event is executed, which is undesirable due to the high communication overhead costs entailed by this operation. However, here we take advantage of the fact that bulk platelets (i.e., free platelets that are not a part of a clot) exist in low concentrations in flowing blood (<1% of occupied LKMC sites). We therefore treat bulk platelets as tracer particles that interact only with aggregated platelets comprising a clot but are otherwise noninteracting with other bulk platelets. In other words, a platelet assumes an excluded volume only after it becomes a part of a clot via a bonding event. This assumption is reasonable (and later will be quantitatively validated by comparison to a serial reference simulation) because the probability is low that any two bulk platelets on different cores execute simultaneous hops resulting in an overlapping configuration. To minimize the impact of overlapping bulk platelets, a periodic scan is performed for platelets occupying the same lattice site, and if found, one of the overlapping platelets is deleted from the system. At any instant of time, the previously known clot configuration is available to all cores. We define an update frequency, F, which is the number of LKMC moves after which a global call is initiated to update the clot configuration in each core. Additionally, we only make global calls to update Γ max after every F LKMC moves by specifying Γ max to a value slightly higher than the highest total rate to account for any fluctuations in the highest total rate (e.g., a platelet entering the system will increase the total rate). The choice of F therefore determines the trade-off between accuracy and computational cost of the algorithm. For all the simulations presented in this paper, we set F to be three times the average number of bulk platelets on a core. In other words, bulk platelets perform on average three lattice hops between successive updates of the aggregated platelet configuration. We maintain good parallel efficiency by making use of these assumptions and approximations to avoid global calls without compromising the accuracy of the parallel model when compared to the results of the serial model (see Section 3).

Soluble platelet agonist concentration fields C _j (x,y,z,t) (where j = ADP, TXA₂) were determined using the FVM solution to the convection-diffusion-reaction equations: ∂ C j ∂ t + v ⋅ ∇ C j = D j ∇ 2 C j + R j (8) where D _j and R _j are the diffusion coefficients and the volumetric release rates of ADP and TXA₂, respectively, and v is the local blood flow velocity. The open-source software package OpenFOAM was used to carry out our FVM computations of spatiotemporal agonist concentration profiles [17]. OpenFOAM is an open-source C++ library that with parallel FVM implementations based on MPI, and has been widely used in both academic and industrial settings. OpenFOAM deploys a spatial domain decomposition approach for parallelization, where the simulation domain is divided among the cores. The mesh and fields were decomposed using the decomposePar utility provided in OpenFOAM. Starting from an initial coarse mesh, the dynamicRefineFvMesh utility provided by OpenFOAM was used to perform topological refinements to the mesh. At mesh locations where platelets become sufficiently activated for dense granule release, the mesh spacing was refined to 1.5 μm. In such a scenario, a single platelet overlapped several FVM cells, so the cell at the center of the activated platelet was treated as the source element for the PDE calculation. The time derivative was approximated using the Crank Nicolson scheme with a time step of 0.01 s.

To model the blood flow velocity profile across the simulation domain, we use the lattice Boltzmann (LB) method, which has been used extensively in blood flow modelling [27–32]. Importantly, LB lends itself well to parallel computing [33]. Specifically, we carry out our LB simulations using Palabos, which is an open-source computational fluid dynamics solver based on the LB method [16]. Palabos has been demonstrated to exhibit excellent parallel efficiencies for benchmark computational fluid dynamics calculations [16]. It is designed in C++ with parallel features using spatial domain decomposition and MPI for communication between neighboring cores.

We assume that blood is an incompressible Newtonian fluid that satisfies the Navier-Stokes equations. However, instead of solving these equations directly, LB models the fluid with fictive particles. The fundamental quantity underpinning LB is the density distribution function, f i x , t , in phase space, x , c i , where t denotes time and c i denotes the lattice velocity along direction i. The evolution of the density distribution function is governed by a discretization of the Boltzmann equation, given by f i x + c i , t + 1 − f i x , t = 1 τ f i e q − f i x , t , (9) where f i e q is the equilibrium distribution function based on the current distribution along direction i and τ is the relaxation parameter. The simplest incompressible Bhatnagar-Gross-Krook (BGK) scheme is used for relaxation to equilibrium via collisions between the molecules of a fluid [34], given by f i e q = w i ρ 1 + 3 c i ⋅ v + 9 2 c i ⋅ v 2 − 3 2 v ⋅ v . (10) The fluid density ρ and velocity v in scaled LB units are computed as ρ = ∑ i f i , (11) ρ v = ∑ i f i c i   . (12)

The fluid kinematic viscosity in scaled LB units, ν L B , is related to the relaxation parameter τ as ν L B = c s 2 τ − 1 2 = 1 3 τ − 1 2 . (13)

In all LB simulations reported in this paper, the fluid domain is discretized using a uniform D3Q19 lattice [ i = 1 , … , 19 ]; see Ref. [33]. The 3D lattice for carrying out our simulations is obtained from a stereolithography (STL) file surface mesh describing the computational domain. Palabos is then used to discretize the domain, i.e., convert the surface description of the domain into a volumetric description by identifying which fluid nodes (or voxels) lie inside the domain. This discretization is done internally in Palabos and is fully parallelized. Moreover, parallel domain decomposition is carried out automatically in Palabos using a homogeneous volume mesh. For all LB simulations presented in this work, the lattice spacing was set to 3 μm. During each LB time step, the density distribution function, f i x , t , is updated according to Eq 9. At locations in the domain where there are bound/aggregated platelets, no-slip (bounce-back) boundary conditions are applied [35]. For more details on implementation of the LB model, the exact form of the equilibrium distribution, and LB unit conversion, see Refs. [16, 33, 34].

The multiscale framework described in this work is composed of several modules that require periodic exchange of information, with individual modules each having their respective open-source libraries or routines. To facilitate the coupling between different modules at each coupling time, we use the Multiscale universal Interface (MUI) [18]. MUI is a C++ library that makes use of non-blocking MPI messages to achieve data exchange between each module with minimal modifications to individual module source codes. The flow of information between modules and a schematic of the time stepping alignment across modules is shown in Figures 3, 4 respectively. LKMC provides the positions of all platelets in the domain and the bonding state of each platelet. LKMC requires the velocity field of the fluid from LB to calculate convective rates of motion and the activation state of each platelet from NN to determine the bonding and detachment rates. The NN provides the activation state of each platelet and the input to the NN requires the concentration fields of soluble agonists from FVM. The LB method provides the velocity field and requires the location of all aggregated platelets from LKMC for the location of the no-slip surfaces. FVM provides the concentration field and requires the release rate of platelets, which depends on the location of platelets obtained from LKMC, activation states from NN, and the velocity field from LB. To interpolate the values of velocity fields at locations required by another module, a trilinear interpolation scheme using field values at the 8 nearest available points is used, while a nearest-neighbor sampler is used for interpolating agonist concentration fields and the mapping of platelet locations to mesh cells in LB/FVM.

Coupling of the individual models occurs after time intervals of Δ t c o u p l i n g = 0.1   s . At the start of the simulation LKMC, LB, FVM, and the NN are all specified by the initial condition for each method. Each method is stepped forward in time until the first coupling time is reached. During this update, all modules share information: LKMC updates the positions of all platelets and the bound states of all platelets in FVM and LB; LB updates the velocity field in LKMC and FVM; FVM updates the concentration field in NN; and NN updates the activation state in FVM and LKMC. This process is repeated until the end of the simulation. The timescale for velocity field relaxation is generally ≪ Δ t c o u p l i n g , i.e., the velocity field reaches steady state in LB significantly before the next update time. Consequently, the LB calculation is performed in a quasi-static mode in which it is evolved only until steady state is achieved for a given particle configuration.

To validate our coupled parallel multiscale model of thrombus growth against serial computations, we considered the dynamics of platelet deposition and aggregation in two case studies. In the first case study we considered a cylindrical domain as shown in Figure 5 (inset). The simulation domain was 0.5 mm long with a diameter of 0.12 mm. A semicylindrical reactive collagen surface patch of diameter 0.12 mm and length 0.25 mm that represents an injury was located at the center of the domain. A constant wall shear rate of 200 s⁻¹ that is typical for venous flow was maintained at the inlet. In the second case study, we considered a 1 mm-long stenotic vessel with an inlet diameter of 0.12 mm. In the central 0.5 mm of the vessel, we introduced a narrowing of the vessel that corresponded to a 75% reduction in flow area as shown in Figure 6 (inset). The stenotic area was assumed to express collagen that extends up to half the circumference and represents the region of injury that triggers platelet aggregation, as shown in Figure 6 (inset). The inlet and outlet of the vessel were maintained at a constant pressure drop that corresponds to an initial inlet wall shear rate of ∼1000s^-1, typical of arterial flow conditions [36]. In both case studies, the bulk platelet concentration was set at 1.5 × 10⁵ μL^-1.

Both case studies were first simulated in serial (i.e., on a single core) and the dynamics of platelet deposition used to validate results from corresponding parallel simulations using 8, 16, and 32 cores. As shown in Figures 5, 6, the platelet counts (defined as the number of aggregated platelets that are part of a clot) as a function of time predicted by the parallel runs agree well with those predicted by the serial run. Also shown in Figures 5, 6 (insets) are snapshots of the platelet aggregates generated in both the serial and parallel cases. Histograms of deposited platelet counts placed in bins along the axial location in the domain showed identical results for simulations using 1, 8, 16, and 32 cores (see Supplementary Figure S2). Consequently, the serial and parallel models can be considered functionally equivalent.

Following validation of our parallel multiscale algorithm, we next considered a benchmark simulation in an idealized cylindrical blood vessel 2 mm long with a diameter of 0.24 mm. We imposed a stenosis with a 75% reduction in lumen area in the central 1 mm of the vessel, as shown in Figure 7A. The inlet and outlet of the stenosis were maintained under a constant pressure drop that correspond to an initial inlet wall shear rate of ∼1,000 s⁻¹. The entire stenotic area was assumed to express collagen and trigger platelet activation and deposition. Again, the platelet inlet concentration was set at 1.5 × 10⁵ μL⁻¹. Under these conditions, we simulated the first 180 s of the formation and subsequent growth of the thrombus using our parallelized code running on 64 CPU cores.

The simulation predicted a rapid rate of thrombus growth after an initial lag time of about 80s from ∼0.1 million platelets to ∼0.5 million platelets at 180 s; see Figures 7A, B. This lag time corresponds to the time at which the first monolayer of platelets adhered to the collagen surface become sufficiently activated to release platelet agonists. The inlet volumetric flow rate dropped to ∼10% of its initial value in 180s due to the flow resistance caused by the growing platelet aggregates. A plot of the aggregated platelet count and the change in inlet flow rate as a function of time is provided in Figure 7B. We also present a snapshot of the deposited platelet mass, the concentration profiles of soluble platelet agonists ADP and thromboxane, and the velocity profile along a central slice of the stenosis after 180⁻s in Figures 7C–F. We note that platelet aggregation is significantly higher near the apex of the stenosis where the observed shear rates experienced by the platelets are highest, a feature that is consistent with prior published computational models and experimental observations [6, 9, 15, 37–39]. This observation is due to the role of VWF in mediating platelet adhesion at high shear rates. At pathological shear rates greater than 3000s^-1, VWF undergoes structural changes from a globular conformation to a stretched conformation, increasing the number of exposed ligands [40–42]. To account for this effect in our platelet adhesion model we included an enhancement factor for the platelet adhesion rates as a function of the local shear rate [9].

The parallel performance of the code was examined using the benchmark case considered in Section 3.2. All simulations were carried out on the Bridges-2 system at the Pittsburgh Supercomputing Center [43]. Parallel performance was analyzed using the observed speedup and the parallel scaling efficiency. For a strong scaling study, the speedup is defined as S K = t 0 t K , where t 0 refers to the simulation time for a reference simulation on N 0 cores, and t K refers to the simulation time for a parallel simulation using K cores. Here, N 0 was defined as 4 cores. The parallel efficiency is obtained from the measured speedup as η K = S K N 0 K .

The strong scaling performance was analyzed at an intermediate simulation time of t = 120s following the initiation of thrombus growth. This time was chosen because it provided a representative snapshot of the system where there is sufficient platelet aggregation leading to a spatially heterogeneous distribution of platelets. We note here that the parallel efficiencies of individual modules and the overall code vary as the simulation progresses in time because they are dependent on the number of platelets considered at any instant for LKMC and NN, the number of LB time steps needed to reach steady flow, and the number of FVM cells. For the system under consideration at t = 120s, the LKMC and NN models accounted for ∼0.2 million platelets, the LB model considered an 82 × 82 × 669 lattice (corresponding to a lattice spacing of 3 µm), and the FVM model considered ∼1.5 million mesh cells. The coupled multiscale code was run to simulate 0.1⁻s of clotting and the time taken by each individual module of the code was recorded. During this period, the coupled solver performed ∼10 million LKMC updates, ∼3,000 LB time steps ( Δ t L B = 1.25 × 10 − 8   s ) until LB achieved steady state, and 10 FVM time steps ( Δ t F V M = 0.01   s ). The observed speedup and parallel efficiency as a function of the number of cores used is shown in Figure 8 for each individual module and the overall coupled code.

The NN module is essentially embarrassingly parallel as it requires negligible communication overhead. Moreover, the parallel particle decomposition LKMC method used to distribute platelets across cores ensures a homogeneous load distribution across different cores. The parallel efficiency of the LKMC module remained at 0.76 despite the spatial heterogeneity in the distribution of platelets. While most platelets are a part of the thrombus and do not appreciably influence the total rate, it is worth noting that they would contribute quite significantly to the spatial heterogeneity in the distribution of rates in the domain. This is because the rates of platelet motion are dependent on the local velocity field which is higher in regions where there is flow narrowing caused by deposited platelets. Consequently, spatial domain decomposition would be a poor choice for parallelizing LKMC. The observed parallel efficiency of the FVM module was 0.53 on 64 cores, while the efficiency of the LB module remained at 0.71 on 64 cores. The parallel efficiency of the overall code was observed to mirror that of the LB module was 0.72 on 64 cores.

A breakdown of the overall computational cost is provided in Figure 9 for the simulation on 64 cores. We observed that the NN and FVM modules contribute negligible computational expense compared to the LB and LKMC modules. The computational expense involved in interpolating and exchanging solutions between different modules using MUI was also found to be negligible. The LB module most significantly dictated the overall computational cost of the model, especially at the later stages of the simulation (at times close to full vessel occlusion) for the present system. This explains why the overall parallel efficiency of the framework was observed to mirror that of the LB module. We note that while these computational times have been obtained for a specific set of physiological and technical parameters, we expect the high parallel efficiency of the overall code to remain robust to changes in platelet density or the mesh resolutions of individual solvers. This is because Palabos (LB) and OpenFOAM (FVM) have demonstrated excellent parallel performance for different benchmark computations [16, 17]. The LKMC code also demonstrated excellent performance over a range of platelet densities in a weak scaling analysis (see Supplementary Figure S4).