Measurements of human blood viscoelasticity and thixotropy under steady and transient shear and constitutive modeling thereof

Human blood is a complex suspension of red blood cells (RBCs), white blood cells, and platelets in an aqueous plasma phase containing dissolved proteins. At rest and low shear rates, the RBCs will aggregate to form linear, coin-stack aggregates called rouleaux, the formation of which is dependent on the presence of plasma proteins such as fibrinogen [1,2]. These RBC aggregates will break up under high shear but will subsequently reform when the shear is reduced, giving rise to a reversible, but time-dependent microstructure [3]. The presence of rouleaux significantly increases the apparent viscosity of blood and gives rise to the complex flow effects that blood exhibits including pseudoplasticity—decreasing viscosity with increasing steady shear; viscoelasticity—both elastic and viscous characteristics; and thixotropy—a time-dependent decrease in the viscosity under shear [4–6]. Due to the weak attractive forces holding the rouleaux together, the yield stress of blood is low (on the order of 1 mPa) and viscoelastic response from the rouleaux network is only observed at low shear rates (<0.1 s⁻¹) [7].

At higher shear rates, the isolated RBCs may also deform under strong flows. This deformation under shear is a current active area of research, and recent computational and experimental works have elucidated interesting shape transitions from a discocyte to a multilobe [8]. Moreover, the motion of the RBC may alternate between a Jeffery orbit and a tank treading motion depending on the RBC properties and flow conditions [9,10]. This complex behavior results in a decreased viscosity from hardened RBCs, a slight shear thinning behavior at moderate to high shear rates (roughly 10–1000 s⁻¹), and viscoelasticity with much faster relaxation times than those associated with the rouleaux [11–14]. Thus, to accurately model physiologically relevant flow conditions for blood, constitutive models that combine thixotropy with viscoelasticity are necessary.

Hemorheology has grown in interest due to both the complex nature of blood and the relevance of blood rheology to disease detection, treatment, and prevention. A number of diseases exhibit discernable changes in hemorheology, including sickle cell anemia, diabetes, hypertension, and heart disease [15–18]. As a result, hemorheological measurements can potentially offer an individualized and straightforward method for detecting the presence of such diseases. However, using hemorheology as a diagnostic requires both data and accurate constitutive models for healthy and disease states. Of equal importance for advancing the engineering design of medical devices is the ability to simulate in silico the in vivo flow of blood, which requires accurate knowledge of the transient material properties of blood. Such simulations have applications ranging from the treatment of atherosclerosis to the development of medical devices [19]. We also note that from a rheology standpoint, blood is one of the most abundant examples of a thixotropic fluid and allows for the evaluation of generalized thixotropy models for a unique and biologically relevant case.

Over the past 100 years, much progress has been made in rheological constitutive modeling of blood. Most early studies only considered the steady shear viscosity of blood and used various generalized Newtonian models with the most popular candidates including the Casson, Cross, Carreau-Yasuda, and Quemada models to represent the flow behavior [20]. While all the generalized Newtonian models for blood rheology incorporate a shear thinning behavior, the most notable difference between the various models is the presence or absence of a yield stress. More recently and partly due to the naturally pulsatile in vivo flow conditions that blood undergoes, increasing attention has been given to the transient rheological response of blood. The earliest and most common attempts to model the transient response are through viscoelastic models. These models are typically modified variations of the Maxwell or Oldroyd-B viscoelastic frameworks [21]. While these models are easy to implement and can describe some of the transient behavior observed at high shear rates, they are insufficient in capturing the transient rheology of blood at low shear rates where the structure is constantly evolving [14,22]. A recent review of many of these models and their ability to describe these data sets were given by Armstrong et al. [22].

To properly capture the bulk response of the evolving microstructure within blood, thixotropy models are necessary. The most common approach to modeling thixotropy is through structural kinetic models. These models rely on a kinetic rate equation for a nondimensional structure parameter that typically varies from 0 to 1, indicating the completely unstructured and fully structured states, respectively [23,24]. However, blood is also viscoelastic and so elasticity needs to be incorporated, where the elasticity is dependent upon this microstructure. One method for introducing elasticity is through distinct elastic and plastic stress components, which are evaluated by splitting the total strain and shear rate into elastic and plastic components [25–28]. Another method is by incorporating the structure parameter with an overarching viscoelastic model, typically taking a Maxwell-like form [29–31]. Such models have been used recently with differing definitions of how the physical parameters, such as the viscosity and elastic modulus, depend on the structure parameter [14,32,33]. For example, extensive efforts have been made to model blood thixotropy by Owens with a model that employs a kinetic rate expression to describe the evolving structure in flowing blood. Owens attempts to relate this expression to the average size of the aggregates, which is subsequently used to describe the bulk response of the sample [34]. We also note that numerous earlier attempts in modeling blood thixotropy were limited by a lack of available experimental data on transient blood rheology, with many of the aforementioned models using data collected by Bureau and coworkers in 1980 [35].

In this paper, we propose a modification to the previously published Horner-Armstrong-Wagner-Beris (HAWB) thixotropy model for human blood shear rheology [14]. The modified model, henceforth referred to as the mHAWB model, reformulates the elastic rouleaux contribution to the bulk flow behavior by unifying the thixotropic structure parameter with a viscoelastic model in a similar approach to the elastic isotropic kinematic hardening (EIKH) model as described by Ewoldt and McKinley [30]. We also show new experimental data for transient step shear change experiments on healthy human blood that greatly expand our previously published hemorheological data set [7,14]. The mHAWB model is fit to a subset of this expanded rheological data set and predictions are tested against additional transient rheological measurements. As we will show, the addition of viscoelasticity arising from the rouleaux microstructure leads to a quantitative ability to capture a nonlinear stress evolution during step-rate change experiments. Moreover, we present amplitude and frequency sweep oscillatory shear data for healthy human blood primarily in the nonlinear viscoelastic region in combination with model predictions to demonstrate the full capability of the model. Lastly, we demonstrate that with the proposed modifications, the model is also able to better predict unidirectional large amplitude oscillatory shear (UD-LAOS) experiments, which are of interest for understanding in vivo blood flow, particularly at low shear rates.

In the following, we present the proposed modifications to the original HAWB in Sec. II. In Sec. III, we outline the experimental methodology and protocol, along with the fitting protocol and methodology for evaluating the model predictions with a focus on the transient step shear change experiments. In Sec. IV, we show the results as well as the new experimental data obtained on healthy human blood and provide a comparison with the original HAWB model for blood rheology. In Sec. V, we provide a discussion of the model and address the present limitations. Lastly, in Sec. VI, we present our conclusions.

As for the original formulation of the HAWB model, the mHAWB model presented in this work assumes the shear stress, $σ12$ is a linear superposition of a structural component, $σ12,R$⁠, and a stress component associated with the isolated RBCs, $σ12,C$⁠,

This formulation follows the assumption of the HAWB model which was motivated by experimental data on nonaggregating RBCs. When no rouleaux are present, the blood sample will still exhibit shear thinning, albeit far less significant than when rouleaux are present [12,14]. Theoretically, this assumption is justified by the need for a limiting viscosity value based on the hematocrit and the presence of a viscous medium. The linear superposition formulation is reminiscent of a Kelvin–Voigt type approach and implicates that the strains and stresses for each of the components are independent. This ensures that viscous effects are always present and feature the appropriate magnitude which is particularly important in limiting cases for the structure parameter. Moreover, it ensures that the aggregation/breakup process occurs independently from the RBC deformation process.

The isolated RBC stress contribution remains unchanged from the original formulation of the model and evolves according to a tensorial extended White–Metzner model [36],

where $γ˙$ is the shear rate and $GC$ is the isolated RBC elastic modulus. The isolated RBC viscosity, $ηC$⁠, is a function of the first normal isolated RBC stress, $σ11,C$⁠, which reduces to the Cross model under steady shear conditions,

where $μ∞,C$ is the infinite shear viscosity, $μ0,C$ is the zero shear viscosity, and $τC$ is a time constant governing the dependence of the viscosity under shear. For a derivation of this part of the model, the reader is referred to the original work [14].

The structural stress component is formulated through a structure kinetics approach where the nondimensional structure parameter, $λ$⁠, evolves according to a kinetic rate equation,

where $τλ$ and $τb$ are time constants governing the overall structure evolution and the rate of structure breakup, respectively. The structure parameter is constrained between 0 and 1 indicating the completely unstructured and structured states, respectively. Note that this formulation differs from the original rate equation for the structure parameter by the absence of a shear aggregation term. Analysis of the HAWB model and data fitting showed this term not to be significant such that model simplification was possible to reduce the number of parameters without affecting accuracy. Brownian aggregation and the shear breakage balance to determine the structure of the rouleaux. It should be clarified that the first term in Eq. (7) is termed a Brownian aggregation term to keep with previous convention, but RBCs in general are typically considered non-Brownian. Instead, this term is more representative of confinement as blood is an extraordinarily dense suspension (typically around 40%). Thus, in developing this model, we assume that the shear driven aggregation effects are negligible as RBCs are always in close contact.

The structural stress is defined as a linear superposition of a yieldlike elastic component, $σe,R$⁠, and a viscoelastic component, $σve,R$⁠,

The elastic stress component depends on an elastic strain which is obtained through a decomposition of the total shear rate into an elastic and plastic component (⁠$γ˙e$ and $γ˙p$⁠, respectively) where the plastic component is defined as

The maximum elastic strain, $γmax$⁠, is linearly dependent on the structure parameter,

where $γ0,R$ is the maximum elastic strain for the fully structured state. Using the plastic rate of strain, the elastic strain may be solved through

which balances the elastic and plastic shear rate with the total shear rate when the structure is not evolving. The differential maximum elastic strain term, included when the maximum elastic strain is decreasing, is incorporated due to the dependence of the maximum elastic strain on the structure parameter. This term assumes that when the structure is breaking up, the elastic strain will follow the maximum elastic strain as necessary. Whereas when the structure is building, the elastic strain will not be changed. Finally, the elastic stress component may be evaluated according to

where the yield stress, $σy$⁠, divided by $γ0,R$ represents an effective modulus for the elastic stress component. This formulation differs from the formulation for the yielding element proposed by Dimitriou and McKinley in that it prevents the elastic rate of strain from ever exceeding the total strain rate—a condition which could be achieved in the Dimitriou and McKinley formulation under certain flow reversal conditions [37].

The remaining stress component, the newly incorporated structural viscoelastic stress $(σve,R)$⁠, is defined according to a Maxwell-like equation with a structure dependent modulus and viscosity,

where $GR$ and $μR$ are the fully structured structural modulus and viscosity, respectively. This representation also assumes that structural breakdown decreases the elastic strain in the material, whereas structural buildup maintains a constant stress. This is incorporated through the additional differential structural parameter term in the case of a decreasing structure. The unification of viscoelasticity with thixotropy used in this work differs slightly from the EIKH model used by Dimitriou and McKinley. In the EIKH, an overall modulus controls the total stress including the yieldlike element as well as viscous element. Subsequently, the elastic strain for the yieldlike element is governed by the viscoelastic plastic shear rate instead of the total shear rate [37]. For the mHAWB model, the viscoelastic element is controlled by the plastic shear rate associated with the yieldlike element. Thus, breaking the yielding network is a prerequisite condition for viscoelastic flow with its own unique modulus. This transition is shown qualitatively in Fig. 1.

Under steady shear, the mHAWB model behaves similar to the original HAWB model with a few notable differences. Particularly, the steady shear solution for the structure parameter

lacks the presence of the shear aggregation term. Nevertheless, the limiting behavior of the structure parameter at zero and infinite shear is preserved. Moreover, the total stress, which arises as a summation of the individual stress components,

features a 1.5 power dependence of the rouleaux viscosity on the structure parameter. A summary of the two models compared in this work, the HAWB and the mHAWB, is presented in Table I.

The handling and measurement protocol employed for the blood rheology results in this work follow the previously established guidelines for blood rheology [7,38]. In compliance with UD’s Institutional Review Board, blood was drawn from the antecubital vein after application of a tourniquet with the donors placed in the seated position. The donors had undergone fasting for 8–10 h prior to withdrawal and had no known health complications. For rheological experiments, 6 ml of blood was drawn into a vacutainer tube containing 1.8 mg/ml ethylenediaminetetraacetic acid. An additional 9 ml of blood was drawn and sent for complete blood count, lipid panel, and fibrinogen activity testing. Selected results of the lab testing are shown in Table II. Both donors fall within typical ranges [39].

All measurements were performed using an ARES-G2 strain control rheometer from TA Instruments equipped with a double wall Couette geometry. The dimensions and measurable range for this geometry may be found in our previous work [7]. The blood was placed in the rheometer 45–60 min from withdrawal, and all tests were conducted within 4 h from withdrawal. Tests were conducted at 37.0 °C which was maintained using a Peltier temperature controller. Shear rates in the device never exceeded 1000 s⁻¹ and a preshear of 300 s⁻¹ for 30 s was used before each test to eliminate memory effects from the previous test.

In this work, experimental results will be shown for steady shear, oscillatory shear, step shear change experiments, and UD-LAOS [28]. For steady shear, the measurements were obtained by holding each shear rate for 40 s or a minimum strain of 50. Between each point, the standard preshear of 300 s⁻¹ for 30 s was applied. As is consistent with previous works for blood rheology, at shear rates at or below 2 s⁻¹, the maximum stress obtained over the transient measurement was recorded as the steady shear value to prevent effects of syneresis—the migration of RBCs away from the walls of the measurement device, which is particularly pronounced at low shear rates [7,38]. At shear rates above 2 s⁻¹, the steady shear value recorded was an average of the last 15 transient stress measurements.

For the step shear change experiments, four conditions were tested in this work: a step to 5 s⁻¹ from 0.1 s⁻¹, a step to 5 s⁻¹ from 300 s⁻¹, a step to 10 s⁻¹ from 0.1 s⁻¹, and a step to 10 s⁻¹ from 300 s⁻¹. To prevent significant syneresis effects for steps originating from 0.1 s⁻¹, the shear at 0.1 s⁻¹ was only held for 30 s—a condition which was also imposed for the model fit. These experiments were conducted using the default hardware filtering setting which automatically selects the best frequency, and no filtering was used for the stress, strain, and shear rate signals. It should be noted that a slight delay period (0.01–0.04 s) was observed for the shear rate. To account for this when modeling, the transient data were shifted, setting the zero time as the time at which the shear rate was within 10% of the final shear rate, and an instantaneous change in the shear rate was then assumed when modeling. This method for shifting the data was determined semiempirically based on the resulting oscillatory shear predictions. The rouleaux time scales present in the model are much longer than the delay and are not significantly affected by this assumption, and although the isolated RBC time scales are shorter, the oscillatory shear predictions confirm that the fit value is appropriate.

For oscillatory shear measurements, both a frequency sweep and an amplitude sweep were conducted. The frequency sweep was conducted at a strain amplitude of 10 and the amplitude sweep was conducted at a frequency of 10 rad/s. For both tests, a delay of 10 s was applied before sampling began. For all oscillation tests, the storage modulus $(G′)$⁠, the loss modulus $(G′′)$⁠, and the relative intensity of the third harmonic $(I3/I1)$ were evaluated. Using the Fourier series representation and retaining only the first and third harmonics for a sinusoidal strain, these properties are defined by the transient stress response according to

where $γ0$ and $ω$ represent the strain amplitude and frequency, respectively [40,41]. The moduli and relative intensity are calculated using the trios software from TA Instruments that analyzes the stress signal through a discrete Fourier transform. The UD-LAOS experiments presented in this work follow the previously established protocol [14] where the shear rate is defined as

In this work, the stress response of these UD-LAOS experiments will be plotted using effective Lissajous-Bowditch plots, where the measured stress is a function of the oscillatory component of the strain defined as the deviation about the average strain due to the steady background flow

A combination of the steady shear flow curve data and the transient step shear change data is used to establish the model parameters, which is then used to predict oscillatory shear and UD-LAOS responses for comparison to additional experimental data. The model contains a total of ten parameters—one fewer than the original model formulation. These parameters are outlined in Table III. Similar to the previous approach to parameter fitting [14], the parameters may be separated into six steady shear (⁠$μ0,C$⁠, $μ∞,C$⁠, $τC$⁠, $σy$⁠, $τb$⁠, and $μR$⁠) and four transient parameters (shaded—$GC$⁠, $GR$⁠, $τλ$⁠, and $γ0,R$⁠), where the steady shear parameters are fit to the steady shear flow curve, and the remaining transient parameters are fit to the transient step shear change data. In agreement with previous works, the fully structured elastic strain, $γ0,R$⁠, was set to unitary as this parameter is best evaluated through startup experiments [28,33]. The error values for the parameters were determined through the standard deviation from 10 parallel optimization runs. The results of these runs are shown in the supplementary material, Sec. S-I [47].

The parameters were fit using a parallel tempering algorithm [42] where the objective function is the sum of the normalized L2 norm of the stress residuals,

The L2 norm is normalized by the number of points in each data set $(Pk)$ and the total number of data sets $(N)$⁠. The resulting optimal parameters fit to the steady shear and transient step shear change data may be found in Table III.

In total, the model contains five ordinary differential equations. These equations may be solved using most common numerical methods. In this work, the equations were solved using MATLAB’s built in differential equation solver ode23 s. The initial conditions were specified as the steady state model conditions at the preshear $γ˙i=300s−1$⁠,

For the oscillatory shear experiments, the quantities $G′$⁠, $G′′$⁠, $G3′$⁠, and $G3′′$ were evaluated through harmonic regression. This can be easily performed in MATLAB by solving the linear equation $Ax=b$⁠, where $A$ is a matrix containing the transient harmonic data (only the first and third harmonics were included here), $b$ is a vector containing the transient stress measurements divided by the strain amplitude, and $x$ is a vector which will contain the moduli for each harmonic.

For brevity, only experimental data from donor 1 will be presented in this section. The experimental data for donor 2 may be found in the supplementary material, Sec. S-II [47] and the results are comparable to those described here. We note that both data sets are also available numerically in the supplementary material [47] in a form convenient for use in model validation. The steady shear data and model fit for both the mHAWB model and the original HAWB model are shown in Fig. 2. Both the mHAWB model and the original HAWB model can fit the steady shear stress data, as expected. The only observable differences between the two models regarding the steady shear behavior is a slightly earlier breakdown of the rouleaux structure with these modifications. This difference arises due to the absence of the shear buildup term in the mHAWB model.

To fit the remaining transient parameters, data from four separate step shear change experiments were used. In these experiments, the shear rate was instantaneously changed to moderate shear rates of 5 or 10 s⁻¹ from a high shear rate of 300 s⁻¹ and a low shear rate of 0.1 s⁻¹. These conditions were selected as the previous shear rates correspond to known conditions (nearly completely unstructured and nearly fully structured, respectively) and the subsequent shear rates are at conditions where a moderate structure is still present in the sample. The results of these tests and the model fits are shown in Fig. 3. For experiments originating from 0.1 s⁻¹, the transient stress profile exhibits an overshoot. This overshoot can be captured by the mHAWB model but is not properly captured by the original HAWB model. This is due to the added viscoelastic term for the rouleaux. In the original HAWB model, the stress resulting from the structure was decomposed into a separate elastic and viscous term. Thus, as blood is shear thinning, any overshoot must be due to the viscoelastic term associated with the RBCs, which is not sufficient to capture the observed overshoot. This experiment clearly illustrates the need for an additional viscoelastic stress associated with the rouleaux, which is now incorporated into the mHAWB model.

For the step shear change experiments originating from 300 s⁻¹, an undershoot is observed at very short times. This undershoot results from the much faster viscoelastic process associated with the isolated RBCs. Subsequently, the stress increases as a result of the structure reforming. It should also be noted that a slight overshoot is then observed, which is not captured by either the mHAWB model or the original HAWB model. This slight overshoot is anticipated to be the result of syneresis, which occurs in blood and exhibits a similar transient stress behavior that is more pronounced at low shear rates [7,38,43].

Using the parameter values determined by fitting the steady shear and step shear change experiments, additional dynamic tests can now be predicted by the model. The model predictions for oscillatory shear data are shown in Fig. 4 and compared against new experimental results. Both a frequency sweep (at a strain amplitude of 10) and an amplitude sweep (at a frequency of 10 rad/s) were conducted. Due to the weak structure within blood, small amplitude oscillatory shear measurements are difficult to accurately obtain. Thus, most of the measurements conducted are in the nonlinear region. Historically, these experiments on human blood have not been well fit by existing transient models [22]. However, the mHAWB model can successfully predict the oscillatory behavior well over the range of conditions. Additionally, the mHAWB model improves over the original HAWB formulation, which predicted both a negative storage modulus at low frequencies and a significantly lower storage modulus at low strain amplitudes.

To further explore the onset of nonlinearity, the relative intensity of the third harmonic is also plotted on the right axes in Fig. 4. For the frequency sweep data shown in Fig. 4(a), $I3/I1$ exhibits a monotonic decrease with frequency. This lower degree of nonlinearity at high frequencies can be explained by the transience associated with the isolated red blood cell component being the main contributor to the bulk response, such that only a weakly nonlinear response is detected. This behavior can be anticipated because with increasing frequency at fixed amplitude the sample is experiencing higher shear rates and less structure is present, as evidenced by the structure parameter curve in Fig. 4(c), so the sample becomes less thixotropic. The modified HAWB model improves the prediction of this nonlinear response.

Unexpectedly, for the amplitude sweep shown in Fig. 4(b), the experimental data show a highly nonmonotonic growth of the relative third harmonic with increasing amplitude, with a pronounced maximum around strain amplitude of one. A significant increase in the relative intensity of the third harmonic arises as the sample transitions out of the linear regime corresponding to the onset of structural breakdown, which is not captured by either model. Following this, a second maximum is observed at much higher strain amplitude, which is overpredicted by the models, but less so for the mHAWB model. These maxima are proposed to correspond to the two distinct, nonlinear phenomena occurring in the blood (i.e., breakdown of the rouleaux structure and RBC deformation). This proposition is supported by the structure parameter curve shown in Fig. 4(d), which demonstrates that the structure parameter is small and nearly constant at the amplitude corresponding to the second peak while it is strongly decreasing at the amplitude corresponding to the first peak. It should be noted that the apparent nonmonotonicity for the HAWB structure parameter curve in Fig. 4(d) is likely a result of the preshear and the failure to reach alternance in the allowed time at the low strain amplitudes. Additional data and model predictions, probing the onset of nonlinearity at differing frequencies, are also presented in the supplementary material, Sec. S-III [47]. Semi-quantitative predictions are observed for the oscillatory third harmonic data for both the models, with slightly better agreement observed with the mHAWB model for the double maxima in the amplitude sweep response and the shoulder in the frequency sweep response. Possible reasons for the inability to capture this significant, first maximum will be addressed in Sec. V.

To further evaluate the model’s ability to capture the oscillatory shear flow behavior of blood, the Lissajous–Bowditch elastic projections for various amplitudes and frequencies were also obtained and predicted as shown in Fig. 5. The corresponding viscous projections may be found in the supplementary material, Sec. S-IV [47]. The original HAWB model predictions for LAOS and UD-LAOS are also included in the supplementary material, Sec. S-V [47]. To be comparable to the data presented in Fig. 4, 10 s were allowed for the curves shown in Fig. 5 to reach alternance. After this point, two cycles or 10 s are shown. The Lissajous-Bowditch projections show a primarily viscous nature with the most significant elasticity being observed at small amplitude and high frequency. Additionally, the most significant nonlinearities arise at low frequency conditions as visualized by the deviation from an ellipsoidal stress representation.

In addition to considering the model’s performance under standard oscillatory shear, the mHAWB model was also used to predict UD-LAOS experimental data as shown in Fig. 6 through the effective Lissajous-Bowditch elastic projections. The UD-LAOS effective Lissajous–Bowditch viscous projections may be found in the supplementary material, Sec. S-IV [47]. Interpretation of UD-LAOS experimental data had been discussed in our previous work [14]. The test offers a unique way to probe the thixotropic and viscoelastic characteristics of a material under large deformations without flow reversal. For UD-LAOS, elasticity will manifest in the effective Lissajous–Bowditch projections as both a vertical shift and a left-right asymmetry, whereas structural breakdown will primarily be observed through top-bottom asymmetry. The mHAWB model can accurately predict the response of the sample to UD-LAOS for a variety of amplitudes and frequencies and improves over shortcomings of the original HAWB model by predicting elasticity in excess of the yield stress [14]. Additionally, the experimental results align with those of standard oscillatory shear in that elasticity is most noticeable at high frequencies and low strain amplitudes, and thixotropy is most evident at low frequencies.

These physically motivated modifications to the HAWB model significantly improve the model’s fit and prediction to the transient experimental data for a range of flow conditions. This was shown through the comparison of step shear change fits and oscillatory shear predictions in Figs. 3 and 4, respectively. To quantify the improvement, the objective function values for steady shear and the step shear change experiments are listed in Table IV for the mHAWB model and the original HAWB model. As noted, both models adequately represent the steady shear data and similar objective functions are obtained. For the step shear change experiments, the mHAWB model decreases the objective function value by 46%—a significant decrease. Additionally, the mHAWB model reduced the number of total parameters from 11 to 10. This overall model improvement can be quantified using the Bayesian information criterion difference $(ΔBIC)$⁠, which is evaluated using the residual sum of squares for the respective models, $RSSi$⁠, as

where n is the total number of data points used in the fitting and $ki$ is the number of parameters in the respective model [44]. The combination of parameter and objective function value reduction results in a $ΔBIC$ of 276 in favor of the mHAWB model for the combined steady shear and step shear change fits.

Qualitatively, the mHAWB model can accurately fit the transient step change experiments presented here that the original HAWB model is unable to. Specifically, in the step shear change experiments from a low shear rate to a moderate shear rate shown in Fig. 3, the mHAWB model can capture the full magnitude of the stress overshoot. This overshoot arises because of the added viscoelasticity associated with the rouleaux. For a purely thixotropic response, a monotonic decrease in stress would be expected. However, the original model predicted a slight overshoot, which arose as a result of the isolated RBC viscoelasticity (occurring over a shorter time scale) and the evolving modulus value for the yield element (admittedly a more phenomenological approach). The requirement for enhanced viscoelasticity associated with the rouleaux can also be highlighted in the UD-LAOS experiments shown in Fig. 6. With minimal viscoelasticity, the HAWB model cannot predict stresses in excess of the yield stress at zero shear rate [14]. Thus, by incorporating rouleaux viscoelasticity the mHAWB model remedies this shortcoming of the original HAWB formulation. For the oscillatory shear data shown in Fig. 4, the most notable improvement is seen at low frequencies. The original HAWB model predicts a negative transient storage modulus, which for thixotropic samples is not theoretically impossible and is not observed experimentally for human blood. The mHAWB model remedies this through the added viscoelastic formulation for the rouleaux.

While the mHAWB model shows promising results, several considerations should be addressed regarding the interpretation of the results. Step-rate experiments show transience due to instrument limitations which require consideration in model fitting as explained in greater detail in the supplementary material, Sec. S-VI [47]. Accurately defining the point of the step change is important in determining the overall structure rebuild time constant and the elastic modulus values. However, the structural evolution process occurs over long enough time scales and should not significantly be affected by our assumptions for the step shear change experiment. The isolated RBC elastic modulus value is more dependent on these assumptions, although the value obtained appears to be sufficient based on the oscillatory shear predictions. In addition, inhomogeneous effects such as syneresis and sedimentation may be affecting the results particularly at low shear rates. Syneresis was accounted for in the steady shear data by using the maximum stress measured as the steady state value at low shear rates as traditionally done in the literature [38]. However, syneresis is difficult to adjust for in other rheometric tests. Furthermore, as blood is a living fluid, it also ages when taken out of the body. These effects in addition to protein adsorption at the interface and thermal history effects were discussed in a recent paper, and it was found that steady shear stress measurements can change roughly 10% in 4 h from withdrawal [7]. However, the changes in the dynamic moduli may be higher. Due to the complexity and potential nonmonotonicity of these changes, these effects were not accounted for in this work.

Intriguingly, double maxima are observed in the third harmonic during the amplitude sweep data shown in Fig. 4(b). For the experimental data, two distinct peaks are observed in the third harmonic relative intensity. These peaks can be associated with the two hypothesized nonlinear processes—the structural breakdown of the rouleaux at lower strain and the stretching of the RBCs at higher strain, where the rouleaux are broken down. The mHAWB model and the original HAWB model also predict two peaks in the third harmonic relative intensity; however, for both cases, the peaks are close together giving the impression that the RBC stretching process dominates. This shortcoming of the model could likely be improved by limiting the influence of the RBC stretching process through a smaller $μ0,C$ value and adjusting the characteristic strain amplitudes that the two nonlinear processes occur at through $τb$⁠, $τC$⁠, and $γ0,R$⁠. However, these restrictions on the model parameters would result in worse fits for the steady shear and step shear change data. Note that there is considerably more noise associated with the measured values of the third harmonic relative intensity and syneresis may be occurring at these flow conditions, which could explain the disagreement. The presence of syneresis is supported by the stress data presented in Fig. 5(b), which shows that the model overpredicts the data, as would be expected for regions in which syneresis is occurring.

Due to the linear superposition of the stress components, the individual contributions to the total stress can be extracted to enable insight into their relative influence. This is shown for the steady shear data in Fig. 7. The stress decomposition for the step shear change experimental fits, LAOS Lissajous–Bowditch predictions, and UD-LAOS Lissajous–Bowditch predictions are also provided in the supplementary material, Sec. S-VII [47]. From the stress decomposition, we can see that the isolated RBC contribution primarily dominates over the range of shear rates. However, at shear rates below 1 s⁻¹, the rouleaux stress contribution becomes most significant, and at shear rates below 1 s⁻¹, the yield stress contribution surpasses the rouleaux stress.

For the oscillatory data, the stress decomposition is shown in Fig. 8. At low frequencies, we can see that the storage modulus tends to be dominated by the rouleaux element and is responsible for the observed plateau that begins to approach the expected behavior of a viscoelastic solid (although we are reminded that this is in the large amplitude region and additional harmonics are relevant). For the amplitude sweep data, again the isolated RBC contribution dominates over most of the measured range. As this component contains one transient parameter $(GC)$⁠, it is likely that $GC$ may be evaluated with more accuracy and independently by oscillatory shear. This is confirmed through a sensitivity analysis on the transient parameters in the supplementary material, Sec. S-VIII [47]. While outside the scope of this work, a more systematic study is required to identify independent tests that can isolate other parameters. At low amplitudes, the viscoelastic rouleaux contribution to the stress becomes more relevant. For the range of oscillatory conditions considered, the yielding element plays a relatively insignificant role. However, the definition for the elastic yielding strain indirectly affects the viscoelastic rouleaux term through Eq. (13). This analysis is useful for identifying what components in the blood most significantly affect the flow behavior for a given flow condition. Clinically, this is can be relevant for drug targeting to treat hyperviscosity syndromes or artificial blood design. However, it should be noted that future work focused on adjusting physical properties of the blood (RBC stiffness, RBC aggregation tendency, etc.) is needed to confirm that the parameters in the model correlate well to their physical counterpart.

Blood is a very complex fluid containing numerous constituents that can influence hemorheology with only limited knowledge of their influence [2,33]. While our model attempts to bridge the gap between the macroscopic flow behavior and the microstructure, it is to an extent phenomenological and only considers the RBC contributions. More research is warranted to isolate the possible rheological contributions of other constituents such as white blood cells, proteins, and platelets. Nevertheless, constitutive modeling is important for hematology research. Previous authors have proposed blood rheology as a diagnostics tool for detecting diseases that feature blood hyperviscosity symptoms [14,20,45,46]. However, to accurately diagnose such diseases, accurate models are needed to isolate confounding contributions. Additionally, constitutive modeling for blood rheology is important for circulatory simulations. Throughout the circulatory system, particularly at bifurcations and contractions, unique strain and shear fields are present which require careful attention to the rheology to properly predict physiologically relevant variables like the wall shear stress. Constitutive modeling offers an accurate and computationally inexpensive way to address this.

New experimental data for step-rate experiments and LAOS and UD-LAOS for two donors reveal that there are multiple contributions to the viscoelasticity at low and high shear rates, as well as thixotropy and apparent yielding at low shear rates in human blood. In this work, modifications to the previously published HAWB model were proposed to account for these contributions. The main modification involved the treatment of the rouleaux stress term—adjusting the term from a summation of a viscous and an elastic component to a summation of a viscoelastic and elastic component based on a more physical description of how the rouleaux evolves during shearing. The number of parameters was also reduced from 11 to 10 by eliminating the shear aggregation term in the structure kinetics equation, which was determined not to be significant over the range of parameters explored. The mHAWB model and the original HAWB model were compared to new experimental data for step shear change and oscillatory shear experiments. Quantitatively, the mHAWB model improved the objective function value and Bayesian information criterion significantly over the original HAWB model. This is due to the ability of the mHAWB model to capture a significant stress overshoot that is observed when stepping the shear rate from a low value to a moderate value. This overshoot is a distinctive feature of a thixo-viscoelastic fluid and arises from the initial stretching of the near fully structured state, followed by a subsequent breakdown. As evidenced by the original HAWB model fit, the overshoot cannot be captured by a primarily thixotropic model. It should also be noted that the time scales associated with viscoelastic stretching and thixotropic breakdown are similar and very short, requiring careful rheological data to detect these transient effects.

The mHAWB model also improved over the HAWB model with respect to the oscillatory shear flow predictions. Prior to this, suitable oscillatory shear predictions have never been achieved for human blood using previously established thixotropy models [22]. The need for an additional memory effect, implemented through the rouleaux viscoelastic term, was further highlighted by the amplitude sweep data at low strain amplitudes where the storage modulus becomes comparable in magnitude to the loss modulus. For the original HAWB model, this was severely underpredicted. Moreover, analysis of the UD-LAOS data confirmed this need as evidenced by stress in excess of the yield stress present at zero shear rate. For the frequency sweep data, the original model predicted a negative storage modulus at low frequencies. The mHAWB model addresses this shortcoming by including a structure dependence for the rouleaux elastic modulus that matches the structure dependence of the viscosity in the viscoelastic equation. The mHAWB model was able to accurately predict the observed plateau that occurs at low frequencies. Our new observations of the third harmonic show a highly nonlinear behavior with multiple maxima that we believe reflects two important nonlinear phenomena, namely, the rouleaux structural evolution and the RBC elongation and deformation. This merits more exploration experimentally. Although the model was unable to capture this rich behavior, we believe that adjustments to the relative intensities of the two processes could lead to a more pronounced double maxima predicted by the model.

There are a wide range of valuable applications for robust constitutive models for hemorheology that can be applied over a broad range of transient conditions. Notably, accurate models are required to properly understand the effects of the individual constituents to engineer devices and protocols to use hemorheology as a diagnostic. Moreover, they can be implemented with reasonable computational effort in flow simulations for the circulatory system as well as to model various flow complications of medical interest, e.g., thrombosis, stents, and arteriosclerosis. Although the mHAWB model is scalar and specific to shear flow, it offers an experimentally tested and thorough yet simplistic description of transient shear blood rheology. Moving beyond the scope of blood rheology, the applications of this model are widespread. Complex colloidal suspensions are widely used in industrial applications and are often subjected to unsteady flows. Although this model has only been tested on human blood, there is potential application for this model to other thixotropic and viscoelastic fluids [24,28].

Funding for this research was provided by the National Science Foundation (NSF) through Award No. CBET 1510837. The authors also acknowledge the Nurse Managed Primary Care Center at the University of Delaware STAR campus for support in drawing the blood samples, Quest Diagnostics for performing the lab tests, and the Environmental Health & Safety department at the University of Delaware for their assistance in transporting the samples. Finally, they also thank Dr. Donna S. Woulfe for her assistance with sample collection and her expertise with the biological interactions within blood.

$A$

Transient harmonic matrix used for harmonic regression

$b$

Intermediate variable

$b$

Transient stress vector divided by strain amplitude used for harmonic regression

$BIC$

Bayesian information criterion

$c$

Intermediate variable

$F$

Function value

$G$

Elastic modulus parameter

$G′$

Storage modulus

$G′′$

Loss modulus

$I$

Harmonic intensity

$k$

Number of model parameters

$N$

Number of data sets used in fitting

$n$

Total number of experimental points

$P$

Number of points in each data set

$RSS$

Residual sum of squares

$t$

Time

$x$

Vector containing the modulus values used for harmonic regression

γ

Strain

$γ˙$

Shear rate

η

Apparent viscosity

λ

Structure parameter

μ

Model viscosity

σ

Stress component

τ

Time constant

ω

Oscillation frequency

$0$

Zero shear value or maximum strain for oscillatory shear

$1$

First harmonic

$2$

L2 norm

$3$

Third harmonic

$11$

First normal tensor component

$12$

First-second shear tensor component

$∞$

Infinite shear value

$a$

Constant governing the rate of shear induced structure aggregation

$b$

Constant governing the rate of shear induced structure breakdown

$C$

Red blood cell contribution

$e$

Elastic component

$G$

Constant governing the relaxation of the elastic modulus

$i$

Initial condition

$k$

Summation variable

$max$

Maximum allowable value

$obj$

Objective function

$p$

Plastic component

$R$

Rouleaux contribution

$SS$

Steady shear

$t$

Transient oscillatory

$v$

Viscous component

$ve$

Viscoelastic component

$y$

Yield

$λ$

Overall structure rebuild constant