Journal of Biomedical Engineering and Biosciences (JBEB)

ISSN: 2564-4998

with the initial number densities  and . The transport equations for resting and activated platelets are convection-diffusion-reaction equations of the form

where

with  and  being reaction constants (Table 2). Red blood cells have a shear-dependent effect on the transport of platelets. One way to simulate it is by adjusting Brownian diffusion fluxes.  denotes the diffusion coefficient of platelets, which is enhanced due to shear rate through the following equation [2]

with the Brownian diffusivity  and the coefficient  [2].

The coagulant concentration  represents the lumped effect of all underlying biochemical reactions in the coagulation cascade [2]. In low shear rate areas, there is a production of coagulant at the wall based on the conditions specified at the boundary. The diffusion-reaction equation for the coagulant concentration is

in which  is a reaction constant (Table 2) and  an effective diffusion coefficient of the form

with , which gets enhanced in areas of low cycle-averaged shear rates [2].

The reactive source term in Eq. 12 shows that the coagulant is produced where bounded platelets concentration is sufficiently higher than its threshold value () and the cycle-averaged shear rate is lower than its threshold value (). The decisive production of coagulant concentration (which is assumed to be identically zero, initially), however, happens at the vessel wall via a Neumann boundary condition. The coagulant flux into the domain through the wall is given by

where  is a reaction constant. If the cycle-averaged wall shear stress  is larger than the threshold value 0.2 Pa and the concentration of bounded platelets at the wall is larger then is zero, while otherwise 

The so-called residence time  of liquid components or platelets in the field, is determined by the transport equation

where  is the self-diffusivity of blood. High residence time marks areas where the platelets spend long time [2].

Table 1: Threshold values in the thrombus formation model.

Table 2: Reaction constants in the thrombus formation model.

3. Rheological Modelling

The rheological model of blood as a shear-thinning liquid with a yield stress determines the extra stress as a function of the rate-of-deformation tensor

with the rate-of-deformation tensor  (in Cartesian coordinates x, y for i or j = 1, 2) defined as

Therefore, the fluid does not flow if the stress magnitude is below the yield stress.  Whereas the arguments  and  originate from the invariants of extra stress tensor

and the rate-of-deformation tensor

In these equations, the shear-rate dependent viscosity and the yield stress appear as liquid material properties, which are modelled as laid out in the following two sub-sections.

3.1. Shear-Rate Dependent Viscosity

Shear rheometric data of human blood was measured by a rheometer with a double cylindrical gap geometry for shear rates ranging from  to . The flow curves were obtained for four different haematocrit values (HCT=30%-60%, with a gradation of 10%) from five different test persons. The shear-thinning, inelastic behaviour of blood qualifies the Carreau model

for representing the flow curves. The four quantities , ,  and  are free model parameters to be determined in fitting the model equation to measurement data. The data are shown in Figure 1. The green curves represent the fits with the Carreau model. For evaluating the viscosity at HCT values intermediate to the experiments, a fit of the model parameters was performed. Thus, all four parameters were available as functions of the HCT value. The pink curves correspond to the viscosities calculated with the parameters from the fit for the corresponding haematocrit, as noted next to the curves.

3.2. Yield Stress

The incorporation of a yield stress in the rheological model is based on the mathematical formulation by Merrill [12] as a piece-wise function of the haematocrit value HCT. The model equation for 30%  HCT  50% is

and for HCT > 50% is

 The equations reveal the yield stress in milli-Pascal with the HCT inserted in percent. The model approximation to the experimental data available in the literature is shown in Figure 2.

4. Problem Formulation

Thrombus formation in the case of TBAD was investigated on the 2D geometry of an abstracted aorta given in Figure 3, with the geometrical parameters in Table 3. For the sake of simplicity, only the ascending aorta, aortic arch without carotid arteries, and the descending part of the thoracic aorta were considered. The dissection is realized by splitting up the descending part of the aorta into a TL and a FL with the width  and , respectively. Both lumens are connected via a proximal and a distal tear. The proximal tear is located at the reference line  and the distal one at . Their width is  and , respectively. The ascending part of the aorta has the width  and length . The width of the aortic arch is equal to the one of the ascending part. The radius of curvature is described by . The arch ends at the angle α. B is the width of the aortic wall.

Table 3: Parameters of the geometry.  = 45°

5. Computational Fluid Dynamics

The model equations were solved with the finite volume method using the open-source software OpenFOAM. The 2D geometry of the flow domain was discretized into 16000 hexahedral cells, with a refining grading towards the walls to correctly resolve the gradients on the wall surface. A grid independency study was performed to ensure that a refinement of the spatial discretization would not change the outcome of the value for the volume fraction of thrombus in the FL. For the numerical schemes, a second-order implicit scheme for the time was chosen, gradient schemes with a second-order central difference and divergence schemes were calculated with a Gaussian first or second-order scheme. The non-orthogonality of wall-normal gradients was taken into account by using an appropriate correction scheme. Laplacian terms are calculated linearly limited by a value of unity. For each simulation, first the periodic solution is reached, and then the thrombus formation equations are included. For the inlet boundary condition of the velocity, a parabolic profile taken from [13] was imposed, with the no-slip condition for walls. Further details about the boundary conditions for the thrombus formation equations may be found in [2] and [9].

The heart rate simulated is 60bpm and, therefore, each cardiac cycle is nominally one second. Note though that the time-scale of thrombus formation in the employed model may not be taken by face value [2], and the times discussed in the sequel are not actual time predictions for thrombus growth.

6. Results

The formation of a thrombus is indicated by the parameter  as introduced in section 2. Complete thrombosis for a cell is reached, if  approaches the value 1. Though the value of unity may not be reached, this parameter shows a clear bimodal distribution of values, which are effectively either zero or larger than 0.98. Note though, that the concentration of bounded platelets affects the flow whenever . A measure of the volume fraction of thrombus is defined through the following formula (normalized different from [11]), to quantify the fraction of thrombus in the false lumen,

The rate of thrombosis, subsequently also referred to as thrombus growth rate is

In order to decrease the necessary computation time, the model for thrombus formation is artificially accelerated by making thrombus growth depending on cycle averaged quantities. This results in two different time scales, a haemo­dynamical one and one for the thrombus growth, neither of which represents the physical time scale. The times provided in the following should only be viewed as a relative quantity between the different simulations.

Figure 4 shows the volume fraction of thrombus in the FL as a function of time. The simulations were performed for three different haematocrit values:  HCT = 35% (turquoise), 45% (yellow) and 55% (red). All simulations were run for 200 cardiac cycles to ensure the stagnation of growth for all cases. However, in all cases this stagnation was reached after 100 cycles at the latest.

For all HCT values, the same general trend can be seen in time, with an increase in the amount of formed thrombus. Furthermore, by increasing the HCT values, less thrombus is formed in the FL. In the early stages of thrombus formation up to approximately 20 cardiac cycles, the difference in relative thrombosis for different HCT values is marginal. As time progresses, the difference of the values becomes pronounced. In all three cases, after 73 cardiac cycles, thrombus formation slows down and reaches an asymptotic value, i.e., thrombus formation effectively stops. The final values for the volume fraction of thrombus, which are anti-proportional to the haematocrit value, are 60%, 57% and 54%. The markers in Figure 4 indicate the points with the maximum growth rate. Thus, it represents a turning point in the growth behaviour, from which on the growth rate steadily slows down. The thrombus growth rate over time for the different cases is provided in Figure 5.

From a qualitative point of view, the growth rates all behave in the same way: an early sharp maximum at the very beginning, followed by a steep decrease in growth and dropping into a valley (local minimum), which is followed by a region of constantly increasing growth rate. After it has reached another peak (local maximum), the growth decreases and is constantly decreasing over a long period of the simulation time, until another characteristic point is reached. At this point, the rate of growth diminishes abruptly and subsequently tends towards zero.

An also quantitatively almost identical behaviour applies to all three cases up to about 20 cardiac cycles. The main differences occur between 20 and 73 cardiac cycles. The higher the haematocrit value, the later the local minimum and second local maximum occur. Thus, a higher haematocrit causes a delay in the characteristic points in the growth behaviour. These points are shown in Figure 5; the values are 39s, 47s and 57s, respectively.

Qualitatively, the behaviour of the growth rate parallels the one of the volume fraction of thrombus: the lower the haematocrit, the higher the growth rate. This trend is displayed in form of the time averaged growth rate in Figure 6.

7. Discussion

Due to the fact that the growth behaviour of all cases behaves qualitatively in the same way, the case for 35% haematocrit shall serve as representative for all cases for the purpose of discussing and explaining the interactions of thrombus growth and haemodynamics. Figure 7 depicts snapshots of the progression of thrombosis by means of the volume fraction of thrombus (colour coding refers to the quantity ) at different characteristic points in time during the simulation. The early stage of thrombus formation is shown after 5 cardiac cycles (Figure 7a). Initially, the thrombus is forming within the wake areas above the proximal and below the distal tear. These regions exhibit low shear rates and high residence time and are predestined for thrombosis. The next characteristic behaviour occurs after 20 cardiac cycles (Figure 7b), as soon as those two wake areas are fully covered with thrombus; the growth rate drops to its intermediate minimum (see Figure 5). From this point on, the difference in the rheological properties becomes apparent. After 25 cycles (Figure 7c), the two thrombi in the wake areas start growing further into the FL. The time span between the valley (local minimum) and the following peak (local maximum) can be attributed to the growth of the two thrombus fronts from top and bottom. The growth rate reaches its maximum at the point when the two fronts merge (Figure 7d). For the given geometry, there is no further growth in previously uncovered areas. Hence, until 73 cycles (Figure 7e) the growth continuous only in the FL at the surface of the existing thrombus. The last characteristic point, at which growth begins to stagnate, indicates a kind of stationary state at which the haemodynamical condition at the wall and thrombus front prevent further grow (Figure 7f).

In order to understand the growth characteristics, the cycle-averaged shear rate, which primarily controls thrombus growth in the employed model, is discussed.

In Figure 8, the cycle-averaged shear rate is depicted for HCT 35% at the same six different simulation times, as previously discussed. The threshold value of the cycle- averaged shear rate is , above which the production of coagulant and bounded platelets ceases, compare Eqs. (12) and (6), respectively.

An according maximum value is chosen for the colour coding in Figure 8, to facilitate identifying the areas in which this threshold value is exceeded. No thrombus formation is observed where the cycle-averaged shear rate exceeds he threshold value. 

The shear rates displayed in Figure 8 explain the observed growth behaviour: due to the low average shear rates, thrombus formation starts in the wake areas, located above the proximal tear and below the distal tear (see Figure 8a-c). Once these two areas are covered with thrombus, the shear rates on the surface of the thrombi are too high, so thrombus growth stops there and continues along the wall opposing the tears, where shear rates are small.  When the thrombus covers about half of the region between the tears, the shear rates are high at all surfaces, such that thrombus growth stops (see Figure 8a-c).

The final state of thrombus formation is depicted for the three different cases in Figure 9. Based on the contour for , this confirms that the qualitatively similar growth behaviour is reflected in a similar thrombosis of the FL for the different HCT values.

The results are in agreement with clinical observations, indicating that higher HCT values are correlated with reduced blood clotting and decreasing clot strengths [14-16].

The current model offers a purely haemodynamic and rheological explanation of this correlation. At the outset, only the FL without thrombus is present. The flow rate distributions are very similar at this point, but the higher value for the haematocrit causes a higher viscosity, in turn higher wall shear stresses and thus the threshold values  and , respectively, are exceeded in more areas. The cycle averaged wall shear stresses at this stage are given in Figure 10. The colour coding facilitates the identification of areas in which  exceeds the threshold value of 0.2 Pa, above which the coagulant production at the wall stops (compare Eq. (14)), which constrains the growth of thrombus. This effect is enhanced with increasing haematocrit and results in a slower, delayed and prolonged thrombosis [14].

Presumably, the dependence of the shear thinning on the HCT value additionally enlarges the cycle averaged shear stresses. The reason is, that the shear-thinning behaviour of blood causes the velocity profile to be flattened in the vicinity of the FL center line, causing steeper velocity gradients near the wall. This entails higher shear rates towards the wall. The shear thinning effect likewise increases with increasing haematocrit value (cf. Figure 1), such that the velocity profile will be more flattened at higher HCT. The effect of the shear thinning behaviour, however, has not been quantified so far.

The current discussion of course only explains why the current, strongly simplified model of thrombus growth shows a haematocrit dependence, which is in line with experimental observations. If these rheological effects actually cause, or at least contribute to, the reduced blood clotting at high HCT values may only be answered in future research.

8. Conclusions

This paper presents the application of a haemodynamic-based thrombus formation model for studying the haematocrit dependency of FL thrombosis in type B aortic dissection. A rheological model accounting for shear-thinning and a yield stress is applied for the blood flow simulations. Shear rheometric data of human blood was measured by a rheometer. The yield stress was modelled based on data from the literature. Our study predicted a higher probability of complete thrombosis for cases with lower haematocrit values. This is due to an increase in wall shear stress, viscosity and shear rate in the FL with the HCT values. The results are in agreement with clinical observations on blood clotting and suggest that the negative influence of HCT on blood clotting might be due to its rheological effects on blood viscosity and shear thinning behaviour.

9. Limitations

Even though the model equations were implemented for the general 3-dimensional case, the presented problem in this work was only solved in 2-D. The applicability for the 3-dimensional case must be ensured, as well.

Acknowledgements

The authors gratefully acknowledge funding by TU Graz within the LEAD project “Mechanics, Modelling, and Simulation of Aortic Dissection.” Blood rheometry was carried out at the Center of Biomedical Research of the Medical University of Vienna.

References