Systems Biology Approach for Personalized Hemostasis Correction

The correction of blood coagulation impairments of a bleeding or thrombotic nature employs standard protocols where the type of drug, its dose and the administration regime are stated. However, for a group of patients, such an approach may be ineffective, and personalized therapy adjustment is needed. Laboratory hemostasis tests are used to control the efficacy of therapy, which is expensive and time-consuming. Computer simulations may become an inexpensive and fast alternative to real blood tests. In this work, we propose a procedure to numerically define the individual hemostasis profile of a patient and estimate the anticoagulant efficacy of low-molecular-weight heparin (LMWH) based on the computer simulation of global hemostasis assays. We enrolled a group of 12 patients receiving LMWH therapy and performed routine coagulation assays (activated partial thromboplastin time and prothrombin time) and global hemostasis assays (thrombodynamics and thrombodynamics-4d) and measured anti-Xa activity, fibrinogen, prothrombin and antithrombin levels, creatinine clearance, lipid profiles and clinical blood counts. Blood samples were acquired 3, 6 and 12 h after LMWH administration. We developed a personalized pharmacokinetic model of LMWH and coupled it with the mechanism-driven blood coagulation model, which described the spatial dynamics of fibrin and thrombin propagation. We found that LMWH clearance was significantly lower in the group with high total cholesterol levels. We generated an individual patient’s hemostasis profile based on the results of routine coagulation assays. We propose a method to simulate the results of global hemostasis assays in the case of an individual response to LMWH therapy, which can potentially help with hemostasis corrections based on the output of global tests.


Introduction
The basis of blood coagulation is a cascade of reactions in which coagulation factors participate [1]. The physiological function of blood clotting is to stop bleeding (hemostasis). It is initiated by tissue factor (TF), a membrane glycoprotein exposed on the subendothelium, which is exposed when the vessel is damaged. Coagulation factors form active complexes on the negatively charged phospholipid surface, providing an increase in the enzymatic activity of free factors by 2-3 orders of magnitude. In the body, the source of such membranes seems to be microvesicles, activated platelets and plasma lipoproteins [2]. The resulting reaction of plasma blood coagulation is the thrombin-induced activation of fibrinogen, which polymerizes to form a fibrin clot [3].
The impairment of blood coagulation can manifest itself not only as a standalone disease, such as deep vein thrombosis. Often, it accompanies oncological [4], immunological [5] and infectious [6] diseases; it can reveal itself during surgery [7,8] or hormonal therapy [9]. The most common way to correct hemostasis is to use standard protocols, in which an anticoagulant drug (in the case of thrombotic complications) is administered at a dose corresponding to the patient's weight. However, this approach does not always lead to hemostasis correction, as it does not take into consideration the individual characteristics of the patient, such as the sensitivity to the drug, different initial states of hemostasis and differences in pharmacokinetics. Thus, in the case of individual therapy adjustment, its efficacy and drug dosing are monitored by means of hemostasis tests, both routine (activated partial thromboplastin time (APTT) and international normalized ratio (INR) [10]) and global (thrombin generation assay, thromboelastogram and thrombodynamics [11]). Global assays allow for the resolution of more details of the picture of the hemostasis state and not only show a general shift towards bleeding or thrombosis but also find the impaired part of the hemostatic system [12]. Monitoring hemostasis with such tests during therapy adjustment allows physicians to figure out whether the current drug is suitable for the current patient and what dosage is needed to normalize hemostasis. Such an approach sufficiently increases the treatment quality, but it also increases its cost and requires a specialist who can properly interpret the results of the test and who will be able to make the correct therapy adjustment based on these results.
One of the possible cost-effective approaches towards personalized medicine is to use computational models of blood coagulation. Such an approach was realized in the works of K.E. Brummel-Ziedins' group, where an in silico model of thrombin generation [13] was used to predict in vitro thrombin generation in different groups: healthy volunteers [14], patients with deep vein thrombosis [15], and patients with acute ischemic stroke or transient ischemic attack [16]. Another group using the same computational model of thrombin generation investigated the effect of prothrombin complex concentrate on the correction of thrombin generation in dilution-induced coagulopathy [17] and the effects of a reduced temperature on thrombin generation [18].
The advantages of such an approach are a reduction in the cost of therapy, a reduction in blood collection and the minimization of the stage of drug selection and dose adjustment.
Heparins are widely used in many pathological states, including cardiovascular diseases, surgery and the post-operative prophylaxis of thromboembolic complications. Lowmolecular-weight heparins (LMWHs) have stable pharmacokinetics, with almost 100% bioavailability. Many specialists believe that constant coagulation monitoring during therapy is not required because clinical doses of LMWHs may be corrected based solely on the patient's body weight [19]. However, the heparin concentration in blood also varied substantially in a healthy population and was only partially correlated with the donor's body weight [20]. A more precise dose correction may be required in many scenarios (e.g., low or high body mass index; critical states; renal failure (creatinine clearance < 30 mL/min); when switching from one anticoagulant to another; and age (children and older adults > 75 years) [19,21,22]. The frequency of bleeding and thrombotic complications during LMWH therapy is 3-5% [23].
In this work, we used a model of blood coagulation that described more than a hundred reactions to predict the results of the global hemostasis test for patients with different coagulation disorders undergoing LMWH therapy.
In order to develop an individual coagulation profile of the patient, we estimated the concentration of the main coagulation factors and inhibitors based on the results of routine hemostasis assays, such as APTT and PT. The individual pharmacokinetic profile of LMWH was based on the developed model described herein, which considers its clearance and distribution based on creatinine clearance and the cholesterol level. After this profile was developed, a comparison of simulations of hemostasis correction and experimental global hemostasis assays was carried out, and we found a good correlation between our simulation-based predictions and the experimentally measured hemostasis state of a patient.

Study Population
We conducted an observational study in 12 patients with IE admitted to the City Clinical Hospital named after V.V. Vinogradov, Moscow, Russia. This study enrolled patients older than 18 years old undergoing LMWH therapy. Enoxaparin sodium 4000 IU (and in 1 case, fondaparinux sodium 2.5 mg) was administered subcutaneously once or twice a day.
Blood samples were collected at three time points: 3 h after the injection of LMWH, 6 h after the injection of LMWH, and before the next injection of LMWH (12 h after the injection of LMWH).
The blood of healthy volunteers (N = 5) was collected for the preparation of normal pooled plasma (NPP).
Written informed consent was given by all participants. All patients were evaluated and treated in accordance with the latest European Society of Cardiology (ESC) clinical guidelines [24]. The study was conducted in accordance with the Helsinki Declaration. The protocol was approved by the Local Ethics Committee City Clinical Hospital named after V.V. Vinogradov. The identification of patients has been made anonymous for the publication of the article.
Thrombin generation assay. The kinetics of thrombin generation in plasma was monitored using the hydrolysis rate of the slow fluorogenic substrate Z-Gly-Gly-Arg-AMC by thrombin formed during clotting, as described in [25]. It was performed in platelet-free plasma (PFP) obtained by serial centrifugation at 1600× g for 15 min and at 10,000× g for 5 min. Plasma was diluted with a mixture of the fluorogenic substrate, phospholipids and tissue factor diluted in BSA vehicle in a ratio of 2:1.
Thrombodynamics-4d. The thrombodynamics-4d assay was performed with a thrombodynamics-4d analyzer and kit (HemaCore LLC, Moscow, Russia) PFP. This method is based on registering simultaneous fibrin spatial clot growth (via a light-scattering signal) and thrombin spatial distribution (via thrombin-induced hydrolyzation of the fluorogenic substrate) after the activation of clotting in a thin layer of plasma after contact with an immobilized tissue-factor-bearing surface [26]. In this experimental setup, lyophilized reagents were used, and thus, blood plasma was not diluted.

Statistical Analysis
We used the Mann-Whitney U test (non-related samples), which is able to estimate the difference between two independent groups with non-normally distributed data in terms of the level of the parameter, measured quantitatively. It allows for the identification of differences in the value of the parameter between groups of small size. Results were considered significant if p < 0.05. Statistical analysis was performed using Origin Pro 2018 (OriginLab Corp., Northampton, MA, USA) software.

Detailed Model of Blood Coagulation
This model describes the biochemical reactions of blood plasma clotting. The chemical design of the model is based on our previous model [27], with modifications. The mathematical model of blood coagulation consisted of 59 volume variables and 19 surface variables for the spatially inhomogeneous case and 78 variables for the homogeneous case. Model parameters, constants and equations are presented in in the Supplementary Materials. An important feature of the developed model of blood coagulation was taking into account the contributions of various types of lipid surfaces to the formation of the system's response. The model considered 3 types of lipid surfaces: (1) lipids initially present in the blood plasma sample; (2) artificial phospholipids in which the coagulation activator tissue factor is dissolved; and (3) artificial phospholipids added to the sample to enhance the response of the system (thrombin generation test and thrombodynamics-4d).

Pharmacokinetic (PK) Model
The LMWH PK model was based on the assumption that after its subcutaneous introduction, the drug first enters the interstitial space, from where it passes into the bloodstream, where its free form is cleared, while its lipid-bound form is protected from clearance.
LMWH, when administered subcutaneously, enters the interstitial space. There is an exchange with blood, causing a gradual increase in LMWH in the blood. At the same time, it is partially removed from the bloodstream by renal clearance and partially associated with lipids, which leads to a slowdown in its excretion. Model equations and constant values are given in the Supplementary Materials. To take into account the individual characteristics of patients, the following amendments were introduced.

Simulation Output Parameters
During the simulations, the time-dependent spatial distribution of the concentration of each reagent was recorded (distribution of thrombin, fibrin, etc.). The spatiotemporal distribution of fibrin was transformed into the kinetics of the propagation of the clot growth front using the following rules. We assume that clotting occurs as soon as the local fibrin concentration is higher than 50% of the initial fibrinogen concentration [28]. The position of the zone where the fibrin concentration is half the initial value of fibrinogen is the coordinate of the clot growth front. Plotting the movement of the clot growth front coordinate over time, the clot size (CS) as a function of time was obtained. The time point at which CS increased from the zero value was considered the lag time (Tlag). The slope of the CS, linearized in the time interval Tlag + 2 min − Tlag + 6 min, was considered the initial velocity of clot growth (Vi). The slope of the CS, linearized in the time interval Tlag + 15 min − Tlag + 25 min, was considered the stationary velocity of clot growth (Vst). The thrombin spatial distribution is presented as a complex-shaped curve with a prominent peak propagating from the activation site. This peak has an almost constant amplitude (A).

Comsol Parameters of Simulation and Numeric Methods
Simulations were performed with Comsol 6.0 (Comsol, Burlington, MA, USA). In our 1D simulations, we used an absolute tolerance factor of 10 −6 (scaled) and a relative tolerance of 0.01. The time-stepping method was generalized alpha (which contains a parameter, called alpha in the literature, to control the degree of damping of high frequencies) with strict steps. Equations were solved using MUltifrontal Massively Parallel sparse direct Solver (MUMPS). We used an area 4 mm long with an equidistant mesh of 400 nodes.
In our homogeneous simulations, we used an absolute tolerance factor of 0.05 (scaled) and a relative tolerance of 0.005.

Detailed Model of Blood Coagulation: Verification and Evaluation
Native blood contains phospholipid surfaces from the following sources: platelets, chylomicrons, and high-, medium-, low-, and very low density lipoproteins. Since we examined a blood plasma sample obtained from whole blood by double centrifugation, there were practically no platelets in this sample, and their contribution could be neglected. Of the other sources of phospholipid surfaces, procoagulant activity was observed mainly only on very low density lipoproteins [29]. An estimate of the phospholipid surface equivalent for very low density lipoproteins gives a value of the order of 200 µM [30]. This value was used to describe the lipids originally present in the sample.
Tissue factor is a transmembrane glycoprotein and must be incorporated into the phospholipid surface in order to function properly. Therefore, the standard reagent used to activate coagulation in coagulation tests is phospholipid vesicles with incorporated recombinant tissue factor. According to the literature data, in such a reagent, there are about 5000 phospholipid molecules per molecule of tissue factor [31]. These data were used to describe artificial lipids that enter the plasma sample upon the activation of coagulation. Artificial phospholipid vesicles composed of phosphatidylcholine and phosphatidylserine (3:1) are widely used in various clinical and scientific studies of hemostasis to enhance the signal and are a necessary reagent for thrombin generation and thrombodynamics-4d tests. Their standard concentration is 4 µM. The main reactions of blood coagulation occur on phospholipid surfaces, where complexes such as external and internal tenase and prothrombinase are formed. In this case, the proteins forming the complex occupy a certain place on the surface, and if this surface is too small, then the formation of complexes will be difficult. In this work, this was taken into account: on each of the three types of lipids present in the model, the formation of complexes of internal tenase (a complex of factors IXa:VIIIa) and prothrombinase (a complex of factors Xa:Va), as well as the binding of prothrombin to the lipid surface, is possible. On phospholipids containing tissue factor, the formation of an external tenase (complex of factors VIIa:TF) is also possible. The binding of factors to the lipid surface reduces the free (reactive) area, thereby slowing down the binding of other factors to lipids.
The model was verified in a homogeneous formulation using the thrombin generation test ( Figure 1). For verification, we compared parameters such as T1/2 (time to reach half of the maximum signal from the fibrin clot), Tmax (time to reach the maximum thrombin concentration) and Amax (maximum thrombin concentration). The thrombin generation test and thrombodynamics-4d test measure the magnitude of the fluorescent signal from a substrate that is cleaved by thrombin [26]. This substrate can be cleaved not only by free thrombin but also by the thrombin-alpha2-macroglobulin complex. Thus, in the test, we registered a signal from two components. In the model, we The thrombin generation test and thrombodynamics-4d test measure the magnitude of the fluorescent signal from a substrate that is cleaved by thrombin [26]. This substrate can be cleaved not only by free thrombin but also by the thrombin-alpha2-macroglobulin complex. Thus, in the test, we registered a signal from two components. In the model, we can separate the contribution of each component, but in the future, to compare the simulation results and experimental data, we will use the total signal of free and alpha2macroglobulin-bound thrombin with a factor of 0.6, because bound thrombin is somewhat worse at cleaving the substrate than free thrombin [32]. Figure 1C shows the integral contribution of thrombin and the thrombin-a2-macroglobulin complex, where thrombin actually means the sum of free thrombin and thrombin in complex with a2-macroglobulin. Therefore, panel D also shows the sum of free thrombin and thrombin in complex with a2-macroglobulin.
The simulations differ from the experimental data in the clotting onset and the signal from the thrombin-alpha2-macroglobulin complex. The first strongly depends on the preanalytical conditions [33], which are not described in simulations. The second depends on the endogenous level of alpha2-macroglobulin and may be different in the experimental and simulation setups.
In a spatially distributed setting, the model was verified using the thrombodynamics-4d test (Figure 2).  It can be seen that the simulation results in the homogeneous case are in qualita agreement with the experimental data, while in the spatially inhomogeneous case, a g quantitative agreement is observed (the parameters differ by no more than 15%). It can be seen that the simulation results in the homogeneous case are in qualitative agreement with the experimental data, while in the spatially inhomogeneous case, a good quantitative agreement is observed (the parameters differ by no more than 15%).
Afterwards, we used experimental data derived from our previous work [34] (Table 1a) to examine the model's ability to describe individual differences in hemostasis associated with differences in the concentration of coagulation factors using data on thrombin generation and thrombodynamics assays performed in factor-deficient plasmas (Table 1b). The mean error for the thrombin generation test (for all parameters and deficiencies combined) is 38 ± 50%, and for thrombodynamics-4d, it is 54 ± 55%. As we can see, the simulation and in vitro measured parameters differ by no more than 1.5 times. The discrepancies in experimental and simulation tests may originate from the uncertainty of coagulation factor concentrations in deficient plasma, as only the deficient factor level was estimated. Another source of the difference may be the amounts of preactivated factors (i.e., cold activation of factor VII [35] and contact activation of factor XII [36]) generated in deficient plasmas during production and storage.

PK Model of LMWH: Development and Evaluation
We found that, according to the kinetics of LMWH clearance, patients can be divided into two groups ( Figure 3): those whose heparin concentration dropped by more than 60% of the maximum (measured at point 1) and those whose heparin concentration fell by less than 60% of the maximum during the same time. deficient plasmas during production and storage.

PK Model of LMWH: Development and Evaluation
We found that, according to the kinetics of LMWH clearance, patients can be divided into two groups ( Figure 3): those whose heparin concentration dropped by more than 60% of the maximum (measured at point 1) and those whose heparin concentration fell by less than 60% of the maximum during the same time.  LMWH clearance in group 1 was significantly higher than in group 2 (Mann-Whitney test, p = 0.05). According to the results of the comparison of these patients, it turned out that, in group 1, total cholesterol was significantly lower than in group 2 ( Figure 3B) (Mann-Whitney test, p = 0.05). The direct comparison of normalized anti-Xa activity and the total cholesterol level revealed that they correlated with r = 0.67 at time point 2 (6 h after LMWH administration) and r = 0.79 at time point 3 (12 h after LMWH administration) ( Figure S1 in the Supplementary Materials).

Pharmacokinetic Model of LMWH
We used individual patients' parameters (Table 2) to develop a PK model of LMWH clearance. LMWH is administered subcutaneously, and the drug first enters the interstitial space, from where it passes into the bloodstream. By analogy with models describing the dynamics of glucose changes [37], several compartments were considered in the LMWH model ( Figure 4).
LMWH is administered subcutaneously, and the drug first enters the interstitial space, from where it passes into the bloodstream. By analogy with models describing the dynamics of glucose changes [37], several compartments were considered in the LMWH model ( Figure 4).  LMWH, when administered subcutaneously, enters the Vint compartment. Through exchange with the Vb compartment, its concentration in the blood begins to increase. At the same time, it is partially removed from the bloodstream by renal clearance and becomes partially associated with lipids, which leads to a slowdown in its excretion.
Model equations and constant values are given in the Supplementary Materials. To take into account the individual characteristics of patients, the following amendments were introduced. Vb is the volume of blood flow, calculated based on the sex and values of the weight and hematocrit of the patient according to the formula: Vint is the interstitial volume, calculated based on the value of the patient's weight according to the formula: Here, W is the patient's weight in kilograms, and HC is the hematocrit (%). The rate of association with Ka lipids depends on the level of cholesterol: Here, Chol is the measured level of total cholesterol (mmole/L). The average half-life of LMWH is 4.5-5 h. In the model, the rate of elimination of LMWH depends on the level of creatinine (which describes the patient's renal clearance status): Here, Creat is the measured creatinine level (µmole/L). The model was validated by comparing the simulation results with the measured values of heparin concentration (Figure 5), and the mean error was quite low.

Personalized Hemostasis Profile
The levels of intrinsic pathway coagulation factors were assessed from the results of the APTT test. The levels of factors of the external pathway were assessed according to the PT test. As shown in [38], the intrinsic clotting factor concentrations (VIII, IX, XI and XII) were inversely proportional to the APTT value, while the PT value was inversely proportional to the extrinsic and common factor concentrations (V, VII and X) [39].
Thus, we assumed that the concentrations of factors Fg, VIII, IX and XI changed in proportion to the change in APTT: F is the initial factor concentration, APTT is the value of the APTT test, and 31.43 is the mean value of the normal range of APTT values of our test system.
The concentrations of factors II, V, VII and X changed in proportion to the change in PT: F is the initial factor concentration, PT is the value of the PT test, and 12.09 is the mean value of the normal range of PT values of our test system.
F0 is the normal average concentration of factor F. The antithrombin concentration was measured directly. The value of very low density lipoproteins changed in proportion to the change in the value of LDL. The concentration of LMWH in the sample was calculated using the pharmacokinetic model.

Personalized Hemostasis Profile
The levels of intrinsic pathway coagulation factors were assessed from the results of the APTT test. The levels of factors of the external pathway were assessed according to the PT test. As shown in [38], the intrinsic clotting factor concentrations (VIII, IX, XI and XII) were inversely proportional to the APTT value, while the PT value was inversely proportional to the extrinsic and common factor concentrations (V, VII and X) [39].
Thus, we assumed that the concentrations of factors Fg, VIII, IX and XI changed in proportion to the change in APTT: F is the initial factor concentration, APTT is the value of the APTT test, and 31.43 is the mean value of the normal range of APTT values of our test system.
The concentrations of factors II, V, VII and X changed in proportion to the change in PT: F is the initial factor concentration, PT is the value of the PT test, and 12.09 is the mean value of the normal range of PT values of our test system. F 0 is the normal average concentration of factor F. The antithrombin concentration was measured directly. The value of very low density lipoproteins changed in proportion to the change in the value of LDL. The concentration of LMWH in the sample was calculated using the pharmacokinetic model.
Evaluation: Experimental data vs. simulation. Figure 6 shows the comparison of the distribution profiles of thrombin and the dependence of clot size on time in the experiment and in the simulation for patient 2 at point 1, corresponding to 3 h after LMWH administration. Table 3a-d show the comparison of the parameters of the thrombodynamics-4d test and the simulation results for patient 2 at time points 1, 2 and 3, corresponding to 3 h, 6 h and 12 h after LMWH administration. The main difference between the in vitro experiment and the simulation is that clot propagation is almost arrested after one hour in the patient's plasma, while in silico, it continues without a velocity drop. As the propagation rate depends on the levels of factors VIII, IX and, to a lesser extent, XI, it may indicate that our method for their concentration estimation may not be very accurate, and further improvement is needed. However, as most of the test parameters were evaluated before 60 min of clotting, the discrepancy in the simulated and measured parameters is reasonably small.

Discussion
One of the variants of using computer modeling in diagnostics was demonstrated in the works of K.E. Brummel-Ziedins' group. They suggested that combinations of variations in the concentrations of clotting factors and inhibitors can lead to abnormal clotting, even if all concentrations are in a normal range separately. They measured the levels of coagulation factors, simulated the generation of thrombin in the model, and compared it with the experimental data [14,[40][41][42]. This approach demonstrated that such combinations can indeed explain bleeding and prothrombotic phenotypes.
In a study by another group [43] using the model in [13], the authors analyzed the quantitative effects of blood plasma dilution on the generation of thrombin in the context of intersubject variability. Using data from the LETS study [44], in which the concentrations of clotting factors were measured in 472 healthy volunteers, the authors simulated thrombin generation in undiluted, 2-, 3-and 5-fold-diluted plasma. Dilution caused a decrease in the thrombin peak height, the area under the curve and the maximum rate of thrombin generation. In another study [18] by this group, the effect of a reduced temperature on thrombin generation was investigated. For this purpose, the authors created a set of random temperature coefficients for each of the model parameters and compared the calculated thrombin generation curves with the range in which they were supposed to lie. The authors showed that hypothermia in the range from 31 to 36 degrees Celsius slowed the generation of thrombin, reduced the maximum rate of thrombin generation and had almost no effect on the thrombin peak height or the area under the curve.    The main difference between the in vitro experiment and the simulation is that clot propagation is almost arrested after one hour in the patient's plasma, while in silico, it continues without a velocity drop. As the propagation rate depends on the levels of factors VIII, IX and, to a lesser extent, XI, it may indicate that our method for their concentration estimation may not be very accurate, and further improvement is needed. However, as most of the test parameters were evaluated before 60 min of clotting, the discrepancy in the simulated and measured parameters is reasonably small.

Discussion
One of the variants of using computer modeling in diagnostics was demonstrated in the works of K.E. Brummel-Ziedins' group. They suggested that combinations of variations in the concentrations of clotting factors and inhibitors can lead to abnormal clotting, even if all concentrations are in a normal range separately. They measured the levels of coagulation factors, simulated the generation of thrombin in the model, and compared it with the experimental data [14,[40][41][42]. This approach demonstrated that such combinations can indeed explain bleeding and prothrombotic phenotypes.
In a study by another group [43] using the model in [13], the authors analyzed the quantitative effects of blood plasma dilution on the generation of thrombin in the context of intersubject variability. Using data from the LETS study [44], in which the concentrations of clotting factors were measured in 472 healthy volunteers, the authors simulated thrombin generation in undiluted, 2-, 3-and 5-fold-diluted plasma. Dilution caused a decrease in the thrombin peak height, the area under the curve and the maximum rate of thrombin generation. In another study [18] by this group, the effect of a reduced temperature on thrombin generation was investigated. For this purpose, the authors created a set of random temperature coefficients for each of the model parameters and compared the calculated thrombin generation curves with the range in which they were supposed to lie. The authors showed that hypothermia in the range from 31 to 36 degrees Celsius slowed the generation of thrombin, reduced the maximum rate of thrombin generation and had almost no effect on the thrombin peak height or the area under the curve.
Existing studies use the thrombin generation test to assess the state of the blood coagulation system, which does not take into account the spatial distribution of clotting reactions or the transport of coagulation factors.
In this work, we propose a systems biology approach for personalized hemostasis correction. This approach is based on the combination of a detailed biochemical model of blood coagulation, which was validated to describe thrombin and fibrin formation in global hemostasis assays, such as the thrombin generation test or thrombodynamics, under normal conditions or in the case of coagulation factor deficiency. Coupled with information on the patient's current clotting factor concentrations, estimated from the results of the routine clotting tests APTT and PT, we were able to describe the output of the global tests, thus indirectly estimating the current patient's hemostasis.
The detailed model of blood coagulation employed was verified to describe blood coagulation in a wide range of initial conditions, such as normal concentrations of coagulation factors and their deficiency. In addition, this model included test-specific parameters, such as the dilution of the plasma sample, amounts of artificial lipids, fluorogenic substrate and tissue factor reagent concentration. An ordinary differential equation model was used to describe the thrombin generation assay, where all reactions occurred in the mixed volume, while a partial differential equation model with a diffusion component was used to describe the thrombodynamics-4d assay, where the tissue-factor-bearing surface activated clotting. The model employed three types of phospholipid surfaces that can be present in the blood coagulation assay: very low density lipoproteins (patient-dependent); artificial phospholipids, which can be added to the plasma sample in the thrombin generation and thrombodynamics-4d assays (assay-dependent); and artificial lipids, which are a component of the tissue factor reagent (assay-dependent). However, this model, as well as an in vitro global test that utilizes blood plasma samples to estimate the coagulation state, lacks information about the condition of blood cells, endothelium and vessel patency and is able to consider only disorders of blood plasma clotting, as well as their correction, but in the other cases, different models and tests should be used.
The global hemostasis thrombodynamics assay employs PFP that is not supplemented with exogenous phospholipids to simulate spatial clot growth [45]. In this test, TF is bound to the surface, and clotting in the bulk of plasma occurs only due to procoagulant microparticles present in the sample. The procoagulant microparticle count in PFP can be as high as 3500-10,000/µL and between 0.35 and 1.0 µm in diameter [46]. Considering their average diameter of 0.8 µm and taking into account that one phospholipid head group is 0.7 nm 2 [47], we can calculate that the phospholipid concentration in PFP is about 4-12 nM. In both of our global hemostasis assays, the thrombin generation test and thrombodynamics-4d, the plasma sample is supplemented with 2-4 µM of artificial phospholipids, which is almost 3 orders of magnitude higher than the phospholipid amount from microparticles. The amount of lipids in the TF reagent is about 8000 higher than TF [48,49], which corresponds to a 40 nM phospholipid concentration in the case of the standard 5 pM TF dose in the thrombin generation test in the case in which no artificial lipids are added. Due to this, we neglected the impact of microparticles on blood coagulation in our simulations.
The patient-dependent coagulation factor level was estimated based on the results of routine hemostasis tests (APTT and PT). This approach let us imitate the results of the global hemostasis assays (thrombin generation and thrombodynamics-4d) for each patient with relatively good accuracy.
Analyzing anti-Xa activity in the group of patients, we found that LMWH clearance was significantly lower in the group with high total cholesterol levels. In order to describe the efficacy of LMWH therapy, we developed a personalized pharmacokinetic model of LMWH clearance. Existing pharmacokinetic models of LMWH are population-based (i.e., do not take into account the individual characteristics of patients) [50]. According to the literature, the main route of elimination of LMWH from the body is renal clearance [51]. However, there is evidence that heparin can bind to apoproteins on very low density lipids [52] and apoprotein E [53]. When developing the pharmacokinetic model of LMWH, we proceeded from the assumption that the clearance of LMWH depends on the level of cholesterol, but this assumption needs additional testing in a larger group of patients in a separate study. Our model described the observed level of LMWH in patients from both groups, with high or low total cholesterol levels, quite well. Embedding this PK model in our detailed biochemical model of blood coagulation let us describe the current patient's state under LMWH therapy. However, the limitations of our model should be mentioned. Our estimations of blood and interstitial volumes are not very accurate and do not take into account extreme values of body weight, the body's tissue composition or impaired kidney function.
We found that the highest deviation of the simulation parameters from the in vitro measured parameters was for Tlag and A (value of the peak of thrombin propagating impulse). However, for each patient, the mean deviation, combined for all parameters, was relatively low.

Conclusions
In silico simulations instead of real hemostasis tests may reduce the financial and time costs of treatment and facilitate therapy adjustments, as it may predict the hemostasis state based on the individual parameters of the patient. Our approach, described for LMWH, may be used for other hemostasis-correcting drugs, such as anti-hemophilia drugs, oral anticoagulants, etc. A possible further step can be the combination of the simulation of plasma hemostasis with more sophisticated models of platelet hemostasis, which may improve the prediction accuracy and widen the outcome parameter spectrum.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/jpm12111903/s1. Figure S1: Correlation of normalized anti-Xa activity and total cholesterol level in the group of patients. (A) anti-Xa activity was measured at time point 2 (6 h after LMWH administration). (B) anti-Xa activity was measured at time point 3 (12 h after LMWH administration). Table S1: Initial values; Table S2: Constants of reactions on the TF surface; Table S3: Reactions on the TF surface; Table S4: Inhibition constants; Table S5: Constants of volume reactions; Table S6: Volume reactions; Table S7: Constants of reactions on the TF-containing lipids; Table S8: Reactions on the TF-containing lipids; Table S9: Constants of reactions on the exogenous lipids; Table S10: Reactions on the exogenous lipids; Table S11: Constants of reactions on the endogenous lipids; Table S12: Reactions on the endogenous lipids; Table S13: Diffusion constants; Table S14: Reactions settings; Table S15: Diffusion settings; Table S16: Boundary settings; Table S17: Boundary reactions settings. These tables contain detailed description of the mechanism-driven model of blood coagulation. Table S18: APTT, PT, anti-Xa activity, antithrombin level, parameters of thrombodynamics-4d 3 h after LMWH administration (point 1); Table S19: APTT, PT, anti-Xa activity, antithrombin level, parameters of thrombodynamics-4d 6 h after LMWH administration (point 2); Table S20: APTT, PT, anti-Xa activity, antithrombin level, parameters of thrombodynamics-4d 12 h after LMWH administration (point 3). These tables contain information about results of all routine coagulation tests performed at all time points. Table S21: Constant values contain information on the parameters of the model of low-molecular-weight heparin pharmacokinetics. These data can be found in the Supplementary Materials.