A comprehensive study on different modelling approaches to predict platelet deposition rates in a perfusion chamber

Thrombus formation is a multiscale phenomenon triggered by platelet deposition over a protrombotic surface (eg. a ruptured atherosclerotic plaque). Despite the medical urgency for computational tools that aid in the early diagnosis of thrombotic events, the integration of computational models of thrombus formation at different scales requires a comprehensive understanding of the role and limitation of each modelling approach. We propose three different modelling approaches to predict platelet deposition. Specifically, we consider measurements of platelet deposition under blood flow conditions in a perfusion chamber for different time periods (3, 5, 10, 20 and 30 minutes) at shear rates of 212 s−1, 1390 s−1 and 1690 s−1. Our modelling approaches are: i) a model based on the mass-transfer boundary layer theory; ii) a machine-learning approach; and iii) a phenomenological model. The results indicate that the three approaches on average have median errors of 21%, 20.7% and 14.2%, respectively. Our study demonstrates the feasibility of using an empirical data set as a proxy for a real-patient scenario in which practitioners have accumulated data on a given number of patients and want to obtain a diagnosis for a new patient about whom they only have the current observation of a certain number of variables.


S1 Supporting Figures
-3: Importance of each feature on the prediction of the platelet deposition using the Random Forest algorithm. Bars depict the relative error in the prediction when the labeled feature is not used. For comparison purposes a bar with the relative error committed when using all features is also included. The confidence intervals represented correspond to the median absolute deviation (MAD) divided by the square root of the number of observations.

S2 Multi-layer mass transfer boundary layer model (MBL)
Multi-layer deposition model In our model, given a certain wall flux of platelets N " 1 , platelets first form a monolayer on the substrate with a deposition rate that depends on the available deposition area, i.e. on the number of platelets P 1 already deposited as follows, is the maximum number of platelets in the monolayer, corresponding to 100% of coverage and d p = 2 · 10 −6 m is the diameter of an adhered platelet [3].
Once platelets form a monolayer, we assume that platelets start to deposit on top of the monolayer at a rate Note that P 2 is the number of platelets deposited above the monolayer independently of the number of layers formed, so that the total number of platelets deposited at a given time is P = P 1 + P 2 To model platelet fluxes, we follow a mass transfer boundary layer approach illustrated in Fig. S2-1. We assume that the substrate to which platelets adhere is regular, and generates a mass transfer boundary layer with a deficit of platelets. The concentration of platelets is a continuous variable that depends on two coordinates. We further assume that the flow is laminar and steady and the mass transfer boundary layer thickness is much thinner than the diameter of the perfusion chamber. Under these hypotheses the velocity gradient, or the shear rate γ, within the mass transfer boundary layer is constant (see Fig. S2-1) and platelets are transported towards the active portion of the wall, located at 0 ≤ x ≤ δ , by flow advection and diffusion. Platelets are deposited on the active portion of the wall following a first order kinetic law, According to this, a free platelet on the wall is transformed into an attached platelet.
The governing transport equation for the concentration of platelets can be written as, While Eq. S2-4 can be analytically solved for this particular set of the boundary conditions [1], [2] , the exact mathematical expressions for the concentration distribution and for the wall mass transfer rates are rather complicated. We therefore use the approximate expression for the surface averaged mass transfer rate when a first order chemical reaction occurs on a wall derived by Pallares and Grau [4], where N " C is the surface averaged mass transfer flux or for an infinite reaction rate, which produces a zero concentration of the reacting species at the wall. In this case the max flux [6], This expression of the surface averaged wall transfer rate is the conventional design relation for mass-transfer wall gauges used to measure the wall shear stress. If the mass flux N " C is expressed using a convection mass transfer coefficient, K C , as N " C = K C C 0 , Equation S2-5 can be rewritten as, Equation S2-7 indicates that the effective convection mass transfer coefficient when a first order chemical reaction occurs on a wall, K e f f , can be understood as the inverse of the sum of the mass transfer resistance produced by the finite rate of the chemical reaction, 1 k , and the resistance by convection in absence of reaction, 1 K C . The validity of Equation S2-7 is discussed in detail in [5]. The differences between Equation S2-7 and the analytical solution considering axial diffusion are less than 3% for γδ 2 D > 3.2 · 10 5 and kδ D > 10 5 . Considering the typically extreme values of the parameters in the perfusion experiments we analyze δ = 2.5 · 10 −2 m, D ≈ 10 −10 m 2 s −1 , γ = 1690s −1 and k ≈ 10 −4 ms −1 the previous quantities are γδ 2 D ≈ 10 10 and kδ D = 2.5 · 10 4 and the error of using Eq. S2-7 instead of the analytical solution is of about 0.5%. Under these conditions the boundary layer thickness is about 10 −5 m, two orders of magnitude smaller than the diameter of the perfusion chamber. Another important assumption in our model is the quasi steady state approximation for the platelet fluxes at the wall. The unsteady mass transfer boundary layer for an infinitely fast reaction, k −→ ∞, was considered analytically by Soliman and [5]. If the concentration of the reactive on the active portion of the wall is suddenly set to zero, the time needed for the wall mass flux to reach the steady state is t ≈ 0.8 δ 2 which corresponds to about four seconds for experiments with δ = 2.5 · 10 −2 m, D ≈ 10 −10 m 2 s −1 and γ = 212s −1 . [4] analyzed the time response of the mass transfer boundary layer for a finite reaction rate. According to the correlation given by these authors the time needed for the wall mass flux rate to reach the steady state is also about four seconds for δ = 2.5 · 10 −2 m, D ≈ 10 −10 m 2 s −1 , γ = 212s −1 and k ≈ 10 −4 ms −1 . This time is much lower than the duration of the experiments, of several minutes, and consequently the quasi steady state approximation for the wall fluxes, adopted to obtain the rates at which platelets are deposited, is reasonable. According to Eq. S2-7, the instantaneous flux of platelets towards the substrate to form the monolayer with an instantaneous length of L = δ 1 − P 1 P ∞ can be written as Correspondingly, the instantaneous flux of platelets towards a layer of length of adhered platelets of L = δ P 1 P ∞ is Introducing Eqs. S2-8 and S2-9 into Eqs. S2-1 and S2-2 leads to and respectively. Eqs. S2-10 and S2-11 constitute a set of ordinary differential equations with initial conditions t = 0, P 1 = 0 and P 2 = 0 with two unknown parameters k 1 and k 2 . Eq. S2-10 can be directly integrated to obtain the time evolution of the deposited platelets on the substrate, while we solve Eq. S2-11 with a fourth-order Runge-Kutta method. The integration of Equation S2-10 leads to where we have arranged the different terms to obtain the non-dimensional governing parameters of the problem. Note that D k 1 δ is the inverse of the Dämkholer number Da = k 1 δ D and D γδ 2 is the inverse of the Péclet number based on the shear rate.  It is apparent that according to Eq. S2-10 and for the set of dimensional parameters considered in this example, the deposition of platelets to form the monolayer depends strongly on the kinetic constant. For k 1 = 10 −4 m · s −1 the time evolutions are very similar to those corresponding to infinitely fast reaction kinetics, k 1 −→ ∞ . In general, the shear rate increases the rate of deposition and this effect becomes more evident for the faster kinetics.

Estimation of the kinetic constants form empirical data
As we mentioned previously, Eqs. S2-10 and S2-11 constitute a set of ordinary differential equations with initial conditions t = 0; P 1 = 0 and P 2 = 0 with two unknown parameters k 1 and k 2 . In our analysis, in order to make predictions of the platelet deposition counts for a given set of test parameters, we need to estimate first the kinetic constants using the training dataset. In order to numerically determine the kinetic constants using the MBL model, we assume that k 1 depends only on the type of substrate used in the experiments. For each set of experiments with a given substrate, we then compute the time evolution of P 1 and P 2 (see Eqs. S2-10 and S2-11). We then perform the calculations for several values of k 1 and k 2 in the ranges 10 −3 ≤ k 1 ≤ 10 −8 m/s. For each pair of values (k 1 , k 2 ), we then compute the absolute difference between the predicted value of the total number of platelets deposited and the corresponding experimental value at a given time. For each different substrate, we select the pair of values (k 1 , k 2 ) that minimizes the absolute difference between the measured and predicted values. Note that this model cannot be directly applied to the cases with stenosis because in this case there is a strong the variation of the shear rate along the substrate.

S3 Calculation of shear rate in experiments with 80% of stenosis
Some experiments were performed using an experimental stenosis of 80%. We wanted to use this variable as a quantitative variable for our models (Random Forest and Phenomenological Model). In those cases we computed the shear rate solving numerically the Navier-Stokes equations in the three dimensional domain that emulate the perfusion chamber. For the numerical solution we used the finite volume CFD code FLUENT 14 (ANSYS Inc., Lebanon, USA). The simulation was computed under steady conditions. As boundary conditions we imposed: a) at inlet a parabolic velocity profile corresponding to the shear rate used in the ex-vivo experiments reported in Table 1, b) at outlet a constant pressure and c) at the wall a non-slip condition. We used the commercial CAD software ANSYS DesignModeler geometry tool (ANSYS Inc., Lebanon, USA) to reconstruct the virtual geometry and the computational mesh of the perfusion chamber with different stenosis. We performed a mesh dependence analysis of the average value of the wall shear rate "WSR" in the stenosis implementing three different structured meshes M1 (275000 cells), M2 (535000 cells) and M3 (750000 cells). The error relative to the WSR calculated with M3 between M2 and M3 was equal to 3%. For all calculations we used M2.