Developing Quasi-Steady Model for Studying Hemostatic Response Using Supercomputer Technologies

Formation of the platelet plug represents a primary response to the vessel wall injury, but may also result in vessel occlusion. The decrease of the local blood flow due to platelet thrombus formation may lead to serious complications, such as ischemic stroke and myocardial infarction. However, mechanisms responsible for regulation of thrombus dynamics are not clear. In order to get a deeper insight into the role of blood flow and platelet interactions in the formation of the primary platelet plug we developed a particle-based model of microvascular thrombosis using quasisteady flow approximation. In order to simulate thrombus dynamics at physiologically relevant timescales of several minutes, we took advantage of the supercomputer technologies. Our in silico analysis revealed the importance of platelet size heterogeneity for describing experimental data on microvascular thrombus formation. Thus, our model represents a useful tool for the supercomputeraided computational analysis of thrombus dynamics in the microvessels on physiologically relevant timescales.


Introduction
Despite many decades of both fundamental and clinical research, complications associated with arterial thrombosis -such as myocardial infarctions and ischemic strokes -still represent a greater cause of mortality and morbidity in developed countries [7]. The reason why some of the vessel wall injuries result in the formation of stable occlusive thrombus, which compromises a local blood flow, while the others not -is one the greatest unresolved issues in the field of hemostasis and thrombosis. Thrombus formation at the site of a vascular injury is a highly complex process, which involves platelet adhesion, aggregation, plasma coagulation and takes place under conditions of the blood flow [6]. An enormous number of mechanical and biochemical interactions occurring within the growing thrombus represent a serious challenge for elucidation of the governing mechanisms and justify the relevance of computational models. During the last decade, various models of thrombus formation have been described, each focused on different aspects of the process [1]. However, a clear understanding of mechanisms responsible for complex dynamics of thrombus observed in vivo is still missing. In order to clarify the role of reversible platelet-platelet interactions and platelet-flow interactions in thrombus dynamics, we have developed a novel particle-based in silico model of microvascular thrombus formation. In this paper, we show that our model reproduces complex dynamics of thrombus observed after laser-induced injury of arterioles in mice [11]. In order to simulate thrombus dynamics on the physiologically relevant timescales, we used quasi-steady flow approximation and took advantage of supercomputer technologies. Our approach involves stochastic modeling. Thus, multiple simulations are required for a statistically valid analysis of thrombus dynamics even for a single combination of model parameters. Such a computational demand can be satisfied by the supercomputer which allows performing multiple simulations in parallel. We performed such computations in several weeks using the resources of the Lomonosov-2 supercomputer.
In this paper, we describe the model in the Section 1. Section 2 is devoted to the model validation and Section 3 embodies evaluation of the model performance. Conclusion summarizes the study and points directions for further work.

Materials and Methods
In order to get an insight into complex dynamics of thrombus observed in vivo, including disruptions of external thrombus layers, we developed a particle-based model, which accurately describes platelet-platelet interactions. To address thrombus dynamics in silico on physiologically relevant timescales of several minutes, we restricted the model complexity by considering a two-dimensional case. Our particle-based model considers platelets as discs with certain radius, which may interact with vessel wall (considered as elastic boundary), other platelets and blood flow. Blood is considered as a Newtonian, incompressible fluid, described with Navier-Stokes and the continuity equations.
The model consists of two main modules. The first one describes the particle (platelet) dynamics in a given flow field, which is considered stationary during the given period of time. The second module is responsible for computation of the velocity and pressure fields in case the platelet dynamics, inferred by the first module, results in the change of the thrombus shape.
Platelet dynamics module. The model takes into account two types of inter-platelet interaction: primary reversible (glycoprotein-Ib-mediated) and secondary, which intensifies with the increase of platelet activation (integrins-mediated). The first type of interaction was described by stochastically associating and dissociating springs, while the second type was represented by deterministic Morse potential. Platelet activation was described as a time-dependent process.
Parameters of platelet interaction models were acquired based on experimental data. We inferred stochastic springs model parameters using experiments on platelet rolling over adhesive surface [3]. Parameters of Morse potential were obtained by fitting experimental data on forces between single activated platelets [8]. The Newtonian equations of motion for each particle were solved numerically with a modified version of the Verlet algorithm [12]. The hydrodynamic force, acting on a freely flowing particle, was described by the Stokes law, while forces acting on platelets within the aggregate were calculated using numerical solution of the continuity Navier-Stokes equations, provided by the CFD module. Platelets were randomly generated at the inlet according to local velocity values and mean platelet concentration profile inferred from experimental data on platelet margination effect [13].
Computational fluid dynamics module. If platelet dynamics calculated using the first module resulted in the change of the aggregate structure (e.g. attachment or detachment of the platelet to or from the aggregate, respectively, or significant displacement of platelets within the aggregate), then CFD module was initialized with new boundary conditions. All platelets, which form the aggregate, were considered as solid impermeable barriers with no-slip boundary condition on their surface (contours). In order to account for fluid motion inside the aggregate, the Developing Quasi-Steady Model for Studying Hemostatic Response Using... effective hydrodynamic radius of these platelets was considered two times smaller than the corresponding value used for the calculation of platelet interactions in the platelet dynamics module. The continuity equation and Navier-Stokes equations were solved numerically using the Open-FOAM simpleFoam solver [9], based on Semi-Implicit Pressure Linked Equations (SIMPLE) algorithm, with parabolic inflow boundary condition on the inlet and zero pressure condition on the outlet.
No-slip boundary conditions were chosen for the vessel walls and platelets within the aggregate. The unstructured mesh was generated using a gmsh generator [4]. When primary calculation is performed, the flow velocities and local pressure values are recalculated based on rescaled boundary conditions corresponding to a constant pressure drop, which better describes the hydrodynamic conditions observed in vivo [2].
The computational domain represented a rectangle with the length of 200 µm and height of 35 µm and mimicked a fragment of the arteriole. An injury was modeled as a row of firmly adhered platelets on the lower side of the vessel wall (Fig. 1). The viscosity equaled 10 −3 P a · s, while maximum flow in the inlet was taken 8.75 mm/s in order to provide surface shear rate of 1000 s −1 , a typical scale for arterioles.

Results
In order to validate the model, we performed the analysis of thrombus dynamics obtained in silico and compared it to the corresponding data derived from in vivo experiments of thrombus formation in the presence of coagulation inhibitor [11]. In order to account for platelet size heterogeneity, the radius of each disk was generated from the Gaussian distribution, which fits experimental data on platelet size distribution [10]. We compared the dynamics of this model with a simple model variant where all platelets have the same radius of 1 micron. Figure 1 shows a typical shape of thrombus and typical dynamics of thrombus height in the model with different platelet sizes. within the thrombus. We analyzed the following parameters of thrombus dynamics (Fig. 1): thrombus surface (S max and S min ) and thrombus heights (h max and h min ) before and after disruption, respectively. We also analyzed the relative size of disruptions, defined as [S max − S min ]/S max . The statistics for these parameters with the number of disruption events is given in the Tab. 1. Our in silico results demonstrate that taking into account platelet size heterogeneity results in thrombus dynamics which better reproduces experimental data. Interestingly, this result indicates that variations in platelet sizes lead to decrease in both absolute and relative size of the disrupting parts (embolus size) of the thrombus and thus seem to increase its stability against the flow.

Evaluation and Analysis of Application Behavior
The described approach was implemented in C++ language using OpenFOAM, gmsh and boost libraries. The experiments were conducted on the Lomonosov-2 supercomputer in compute partition (based on Intel Xeon E5-2697 v3 processors with 14 cores). For each parameter set the computations were performed on one node. In order to get the statistics, 14 simulations with different random seeds were initialized, each at a single thread mode, therefore utilizing all available physical cores in the node.
To estimate the increase of calculation time associated with quasi-steady hydrodynamics computations, we compared the mean time required to calculate 60 seconds of physiological time in a stationary model (when the flow is stationary and no CFD simulations are performed) and quasi-steady model. The obtained values for n=10 runs are the following: for stationary model, the calculation time (M ± SE) is 560 ± 50 min.; the time for quasi-steady model is 1960 ± 180 min.
In order to find out the most computationally intensive parts of our implementation (that are mostly responsible for the significant increase in execution time), we have performed program profiling using Intel VTune Amplifier [5]. It should be noted that in these experiments we calculated a much smaller time period -only 500 milliseconds of physiological time. This was done due to two reasons: 1) the profiling process requires a large number of experiments and therefore is very time-consuming; 2) we have chosen an interval of physiological time that reflects a representative as well as computationally expensive part of the simulation.
The results have expectedly shown that the CFD module tends to consume most of the execution time (∼60% in general), leaving only 40% to the platelet dynamics module. The most computationally intensive function in CFD module is responsible for infering the indexes of the edges which surround the given point in space for further calculation of the interpolated values of pressure and velocity -it takes up to 9.4% of the overall execution time. Next time-consuming steps in CFD part are marking particles which belong to the aggregate as well as calculation of distance between particles, since they consume ∼4% and 3.5%, correspondingly.
The interesting part was to discover that a significant part of the execution time is spent on performing useful arithmetic operations -sin/cos, simple arithmetic and exponent calculations (15.1%, 8.7% and 5.1% of execution time, correspondingly). Moreover, performing sin/cos operations is the overall most computationally intensive function. All these operations are called from one major function that calculates the interactions between particles within the aggregate, which makes this function one of the primary candidates for further optimization. Another way of optimization, which seems to be more preferred in our case, is to try parallelizing this program. According to the analysis results, it seems that the most suitable opportunity for parallelization is to perform a recalculation of particle dynamics in parallel (e.g. using OpenMP), since this part is performed in a nested loop with almost independent iterations where the program spends 40% of execution time.

Conclusions
In this paper we describe a novel in silico model of thrombus formation and its implementation, which allows studying thrombus dynamics on the physiologically relevant timescales of several minutes. Using this new instrument we demonstrate the importance of platelet size heterogeneity for describing experimental data on thrombus dynamics and embolization. To perform mass experiments based on this model, we have used computing resources of the Lomonosov-2 supercomputer. In order to find out the performance bottlenecks of our implementation, we have performed program profiling using Intel Vtune software. The analysis results showed us the most computationally intensive parts of our program and allowed us to identify possible directions for optimization.