Recent Advances in Biofluid Mechanics and Bio- and Hemorheology Collating Recent Advances in Predicting Complex Behavior of Human Blood With Thixo-Elasto-Visco-Plastic Models and Sequence of Physical Process

For years industrial polymer production has driven the development of rheological models to characterize the flow of materials. With the evolution of these models has come a corresponding advancement in the understanding of the complex mechanical properties. Recent efforts have been focused on modeling the behavior of complex fluids such as blood, whose microstructure leads to has simultaneous characteristics such as: thixotropy; elasticity; plasticity; and an evolving viscosity (part of which originates with the rouleaux’s evolution). The specific complex behavior of human blood can be analyzed via the analysis of Large-Amplitude-Oscillatory-Shear (LAOS) and Small-Amplitude-Oscillatory-Shear (SAOS) response tests. Unique features of human blood cannot be replicated in legacy steady-state models and, thus, have required the development of more comprehensive models capable of accurately fitting both steady state, transient flow and oscillatory shear flow. Expanding upon prior transient models, collaboration between the Chemical Engineering departments of the United States Military Academy and the University of Delaware has produced a new model, tensorial enhanced structural stress thixotropic-viscoelastic model (t-ESSTV). This model can capture the timescales contained within the plasma and individual red blood cells viscoelasticity and the thixotropic timescales associated with rouleaux breakdown and aggregation. The efficacy of t-ESSTV is demonstrated with a single Donor before consolidating the best fit model parameters of twelve Donor sets of rheological data. We then show the parametric correlations between model and physiological parameters and with the models’ prediction of microstructure, we correlate microstructure with the “elastic, solid-like” metrics as computed by Sequence of Physical Processes (SPP).


INTRODUCTION
As a complex suspension, the aqueous plasma medium of human blood contains three primary structures: platelets, white blood cells, and red blood cells (RBCs). Unique among these structures is the tendency of suspended RBCs to agglomerate into "rouleaux" structures, where the RBCs stack themselves faceto-face. While several factors can affect the formation of these rouleaux, their formation can be dependent upon the concentration of binding fibrinogen proteins and the stress and shear rate applied [1,2]. The binding of the RBCs to one another is reversible, allowing for the repeated break-down and reformation of rouleaux under evolving flow conditions [3]. The structural change of these agglomerates entails a cascading change in blood's general fluid properties. For instance, the viscosity of blood increases upon the formation of new rouleaux due to increased resistance to flow from increased surface area.
In addition to changes in viscosity, rouleaux formation and disintegration affect the thixo-elastic-visco-plastic (TEVP) properties of human blood. Being a complex fluid, containing a variety of suspended rouleaux (microstructure), blood displays the following unique properties; thixotropy, viscoelasticity, and viscoplasticity [4][5][6][7]. Thixotropy describes the evolution of blood's microstructure (rouleaux); viscoelasticity refers to the demonstration of simultaneous viscous and elastic characteristics; and viscoplasticity relates to the presence of a non-zero yield stress threshold above which particles begin to experience irreversible deformation with no elastic structure recovery. The thixotropic and viscoplastic tendencies of human blood can both be primarily attributed to the creation and breakage of the rouleaux structures, which is a function of local concentration of RBCs (or hematocrit). In addition to governing the formation of microstructure, the elasticity of individual suspended RBCs independently contributes to the viscoelastic properties of blood. Nonetheless, the actual yield stress of blood is low at around an average value of 4 mPa maximum, due to the weak intermolecular forces present within the rouleaux microstructure [4,[8][9][10]. This is corroborated by Merrill et al. [4], and Armstrong and Tussing [10].
Of note are the different conditions in which the viscoelastic contributions from either the rouleaux proper or individual RBCs predominate; in higher shear regimes, blood viscoelasticity is primarily contingent upon the elastic deformation of individual RBCs, whereas at lower shear rates the rouleaux play a more prominent role [10,11]. When juxtaposed with the viscoelastic relaxation time contribution of rouleaux, the individual RBC viscoelasticity possess a shorter relaxation time with shear thinning observable shear rates along a range of 10 s −1 to 1,000 s −1 [8,[11][12][13][14][15][16]. This complex dynamic behavior necessitates the development of novel model features able to characterize all aspects of blood's TEVP phenomena to better understand and address blood pathologies [17][18][19][20][21].
While current models can more effectively model the transient thixotropic and viscoelastic behavior of complex flow, firstgeneration generalized Newtonian frameworks such as the Casson and Herschel-Bulkley models allowed for the representation of a spectrum of steady-state properties, only, including shear thinning and yield stress [22,23]. While these models are unable to predict transient and viscoelastic phenomena produced by the evolution of blood rouleaux, they provided the groundwork by which more sophisticated models could be developed, able to represent the microstructure hysteresis and accompanying effects on material elasticity and yield stress [8,9]. The next generation of rheology models, notably exemplified by the Maxwell and the Oldroyd-8 family of models, phenomenologically model viscoelastic fluid behaviors, allowing for the prediction of transient behavior at high shear [24][25][26]. However, to effectively model the behaviors of dynamic, microstructure-laden fluids at low shear, the integration of a thixotropic component is necessary [24][25][26][27][28][29][30][31]. In representing the time-dependent construction and destruction of microstructure, a kinetic rate equation is introduced with a non-dimensional structure parameter varying between zero, representing no agglomeration, and one, representing full microstructure, or fully structured rouleaux [27][28][29][30].
However, representing the microstructure contribution to viscoelasticity requires the simultaneous incorporation dual elastic and plastic stress terms or, alternatively, the combination of the thixotropic component with an established viscoelastic model [31,32].
Total strain and its derivative by time, total strain rate can be separated into two components, allowing for the definition of elasticity and plasticity terms. This subsequently allows for the representation of isotropic hardening (IH) and kinematic hardening (KH) per the kinematic hardening theories of plasticity [28,31,33]. IH describes the relationship between a material's thixotropic behavior and its structural parameters, whereas KH enumerates the effective yield stress dependency upon deformation. The latter of which can also act to retard the back stress response given shear stress [33][34][35][36][37].
While prior models were developed featuring kinematic hardening [38], it was only with the collection of substantial amounts of substantial transient experimental data, rigorously collected and assessed, that new models could be evaluated [8]. In addition to being probed for steady-state characteristics, the comprehensive dataset included response drawn from Unidirectional Large Amplitude Oscillatory Shear (UD-LAOS) experiments, which subjected the blood to both steady and a linearly superimposed oscillatory shear. For UD-LAOS, the minimum possible shear is zero, in contrast to conventional LAOS which can venture into negative shear regimes. The newfound transient data [39], and model fitting allowed for a far more accurate representation of RBC viscoelasticity via the application of a generalized White-Metzner-Cross model, producing the Horner-Armstrong-Wagner-Beris or, as it is now known, HAWB [8]. Whereas this framework itself evolved from the Apostolidis-Armstrong-Beris (AAB) and the Modified Delaware Thixotropic Model (MDTM) [8]. A subsequent version of the HAWB system, the modified HAWB model (mHAWB), integrated the rouleaux viscoelastic response. Other variants, such as ethixo-mHAWB [10], have been presented with slight modifications to the kinetic equation to represent the structural parameter more accurately in the use of a shear aggregation term. A further modification, relevant to this work, manifested in the enhanced thixotropic viscoelastic (ETV) model [9,40].
Related work has explored distinct aspects of various types of rheological models. Instead of the White-Metzner-Cross model as applied to the mHAWB model, other efforts attempted to explore potential advantages to applying modified Oldroyd-8 and Giesekus models to describing the transient and steady-state aspects of complex fluids, including human blood [10,24,41,42]. Armstrong et al. [43] addressed the characterization of blood rheology in modifying the viscoplastic Herschel-Bulkley model while ignoring the viscoelasticity contributions of individual RBCs, as proposed by Saramito [31]. Per Saramito's approached, the entailing ethixo elastoviscoplastic (ethixo EVP) introduced a thixotropic structural parameter bound to a kinetic equation, allowing for a full TEVP representation of blood [43]. Ethixo EVP and other thixotropic EVP variants were found to fit to transient experimental blood data comparatively better than base EVP. Other models included the Multimode-Lambda + Isotropic/Kinematic Hardening or ML-IKH, developed by Wei et al., which combined isotropic kinematic hardening with independent thixotropic time scales to achieve similar modeling capabilities for complex rheology [36]. A similar, modified variant of the Varchanis et al. model, the Saramito-Phan-Thien-Tanner + Isotropic/Kinematic Hardening or SPTT-IKH, combined several elements of prior models to better approximate TEVP behavior [31,37,44]. Later efforts tensorialized these models, possessing anywhere from 11 to 15 parameters.
The developments were accompanied by the parallel development of the viscoelastic enhanced Modified Delaware Thixotropic (VE-MDTM) in incorporating a viscoelastic evolution timescale into the total stress term produced by rouleaux module of the original MDTM system. In this, the presence of a unique timescale associated with the viscoelastic stress contribution of suspended rouleaux was successfully demonstrated [45][46][47][48][49]. Per the Armstrong-Frederick kinematic hardening rule, KH and IH approaches can be applied to currentgeneration TEVP models, allowing for enhanced model step-test performance when used in conjunction with other thixotropic and viscoelastic terms [36,37,43,44]. The adoption of separate timescales for the RBC viscoelastic, rouleaux viscoelastic, and rouleaux thixotropic contributions, featuring additional sheardependent structure destruction and agglomeration timescales, has been found to improve model performance in predicting both steady-state and transient flow. The ethixo-mHAWB, as well as other contemporary models such as that recently introduced by Pincot et. al. [9,40,45], utilize these independent timescales to maximize model performance while still aiming to limit the number of parameters required for operation. These timescales offer the necessary insight into in-vivo rouleaux evolution to permit the future development of means by which to apply advanced CFD modeling to the analysis of blood flow.
The ETV, and ESSTV [46] performance in fitting steady state and transient data to observed blood rheology was substantially improved compared to earlier-generation models; however, the most recent modifications of ETV have applied a novel viscoelastic module to apply recently published approaches to plasticity to better express both the plastic and elastic contributions to stress from the rouleaux microstructure [31,36,37,[46][47][48][49]. In understanding the rheology, and mechanical properties of blood, it is useful to analyze its response under a few experimental shear conditions due to the dynamic nature of blood viscoelasticity. In this, LAOS is a useful analytical tool in investigating its nonlinear behavior through variation of temporal and spatial frames simultaneously. Oftentimes supplementing oscillatory shear analysis is a set of moduli: where γ 0 represents the strain amplitude, ω the oscillation frequency, G′ the storage modulus, and G″ the loss modulus. The storage and loss moduli are appropriate to reconstruct a linear oscillatory shear signal typically found in the small amplitude oscillatory or SAOS regime. Should there be a requirement for higher harmonics to describe the signal, the LAOS regime has been reached. Fourier transforms and the Fourier-Chebyshev derivation approach to elastic and viscous stress are also applied [40,41,43]. Comprehensively describing the derived waveforms allows for their translation to Bowditch-Lissajous plots and Pipkin diagrams. While the Fourier approach has permitted the characterization of the more dynamic aspects of material response, it lacks a complete physical interpretation [41]. However, this can be remedied in the formation of an orthonormal foundation via the use of the first Chebyshev polynomials [42,43]. Recent developments by Rogers et al. have led to the development of the sequence of physical processes framework (SPP), which couples with LAOS to analyze stress response. SPP's titular physical processes relates to the analyzed material's plastic, elastic, viscosity, and yielding properties. Unlike the Fourier and Chebyshev approaches, SPP emplaces unique viscous and elastic moduli along a spectrum of datapoints to calculate transient moduli from LAOS data [50][51][52][53]. TEVP models such as ESSTV can then be functionalized in predicting blood structure through LAOS wit subsequent comparison of those drawn from SPP.
Our main (new) goals here are: 1. to strengthen previous rheological model and physiological model parameters; 2. Show the aggregated rheological model parameters over twelve donors; 3. Show how sequence of physical processes and our predictions of rouleaux are correlated for Donor5, while comparing to the aggregated SPP for the 12 Donors. The remainder of the paper is organized as follows: Section 2. is a brief background, then description of the tensorial ESSTV model. Section 3. Materials and experimental methods, broken into three subsections: materials and data acquisition; parameter optimization; Sequence of Physical Processes description and parametric correlations. Section 4. Experimental and model fitting results. Section 5. Discussion and Section 6. Conclusion and future work. component from the plasma and suspended, deformable, individual red blood cells, shown here in component form: where, Gc is the elastic modulus of the plasma and individual RBCs and μ C (σ ve,xx ) is mathematically described below in Table 1 [8,9,46]. The equations have been derived previously, using a White-Metzner approach, with a Cross model substituting for the viscosity: μ( _ γ) , where μ 0,C is the zero shear viscosity, μ 0,∞ is the infinite shear viscosity and τ C is the Cross constant. We model the rouleaux as 'microstructure' taking values of [0 1], where 0 represents no rouleaux, and 1 represents maximum rouleaux formation at rest ( Table 2). The following structural evolution equation has three terms: 1. Shear induced rouleaux breakdown; 2. Shear rouleaux aggregation; and 3. Aggregation due to Brownian motion.

Dλ Dt
where τ break is the characteristic time of shear breakage, τ aggr is the characteristic time of shear aggregation, and d is set to ½ for human blood, an optimum value from the literature [8-10, 40, 41, 43] (Previous work has incorporated the characteristic time of shear aggregation outside the absolute value; we note this is for convenience and does not change the mathematics of structure evolution [41,46]). We note here that _ γ p is the plastic shear rate and the portion of the total shear rate that is associated with plastic deformation of the microstructure derived from: γ γ p + γ e ↔ _ γ _ γ p + _ γ e [8, 9, 33, 33-35, 41, 45, 46]. Whereby the strain and shear rate are broken into elastic and plastic components. The auxiliary equations of this formalism are listed in Table 1 below, with full descriptions in previous work [8, 9, 33, 33-35, 41, 45, 46].
The demonstrated behavior of rouleaux under evolving force conditions, are the descriptions of rouleaux breakdown and agglomeration over time per τ break and τ aggr respectively [46]. These each present the breakdown and agglomeration of rouleaux in a comparative ratio to the Brownian microstructure aggregation constant. Note that both τ break and τ aggr have units of time, therefore they are used here to define dimensionless Weissenberg numbers controlling the stressinduced rouleaux breakage and stress-induced aggregation speeds as: Wi break τ break _ γ p and Wi aggr τ aggr _ γ p . These are then in turn used to define practical times for stress-induced breakage and reformation as:τ break τ λ /Wi break andτ aggr τ λ /(Wi aggr ) d [46]. One can compare these values at respective shear rates to determine the overall rate of breakage and reformation of the structure. The parameters τ break and τ aggr  [9,41,43,46,51].

Model ESSTV
PAR. μ 0,c , μ ∞,c , τ c , σ y0 , τ break τ aggr , γ 0,R , d, m, G c , G r , τ λ SS σ ss σ ve,ss + μ r λ m ss _ γ + σ y,0 λ ss sign( _ γ)  are incorporated because they have specific physical interpretation and because they both have similar units can be used a comparator for the overall effects on the rouleaux at any given shear rate [46]. The effects of the microstructure will factor into an elastovisco-plastic (evp.) contribution to total stress as follows G R is the rouleaux elastic modulus, μ R is the rouleaux viscosity, σ y,0 is the yield stress and m are a fitting parameter set to 3/2 for blood as per the literature [9,46]. The steady state constitutive equations for the structure and total stress are The remaining auxiliary equations for the t-ESSTV model are shown in Table 1. Per procedures described in the recent work of Horner et al. and Armstrong et al. [8,9,46]. The models can be cast into tensorial form, further enhancing the general predictive performance [46]. The tensorial version of Eq. 6, for the yx, xx, yy, and zz components is shown in Table 3. Please refer to the publications for a detailed description of development procedures [7-10, 40, 41, 43, 46]. In Table 3 below we take the meaning ofσ as follows (from Armstrong et al. [46] which followed the work of Wei et al. [36] whereσ is the equivalent shear stress, "D" is the deviatoric part of the tensor, I is the unit tensor, ":" is the double dot product, and the trace is defined thusly, tr(σ R ) (σ R,xx + σ R,yy + σ R,zz ) [46]. Additionally, we define γ (1) as follows The t-ESSTV model is used here as framework to be applied to the mechanical and rheological interrogation of ten human blood donors. From the best fit parameter values of the t-ESSTV model and physiological parameter values measured in the laboratory, a sequence of parametric correlations is obtained and analyzed. Additionally, we correlate our model predictions for microstructure, λ, with metrics of "solid-like" and "liquid-like" behavior to see if there is a relationship and/or correlation between these disparate measures [46].

Materials and Data Acquisition
Collection and rheological blood analysis was executed via the procedures and protocols established by Horner et al. [8,9,39]. This entailed the drawing of 6 ml of blood into a vacutainer tube with 1.8 mg/ml EDTA. 9 ml of blood was also preserved for subsequent testing via the Complete Blood Count, Lipid Pane and Fibrinogen Activity test, as shown below in Table 4 [8,9,39]. All testing was done under the purview of University of Delaware's IRB.
Data collection was performed with a TA Instruments ARES G2 rheometer, featuring strain control and a double wall couette [8,9,39]. Tests were conducted at 37°C and limited to shear rates below 1,000 s −1 . A preshear of 700 s −1 was applied to negate the effects of thixotropic hysteresis.

Model
Tensorial ESSTV where y i represents a steady state stress datapoint and f i is the relevant model prediction. N represents the number of analyzed points. The second step involves the fitting of transient parameters to three iterations of step-up and step-down shear rate testing experiments. This again calls upon the parallel tempering algorithm: Where M represents the number of step-up or step-down in shear rate tests. Fitting the model to collected rheological data is accomplished via the use of Matlab's ode23 function to process the component differentials [8,9,41,43]. This is done to pinpoint a good initial guess for the parameters during the tensorial model parameter fitting. During the t-ESSTV parameter fitting, using the initial guesses generated above the total F cost value is computed by normalizing both the steady state and transient data with individual stress measurements for each steady state point for the steady state data, and maximum stress values for the entire step-up/down for the step-up/down in shear rate data sets as shown below in Eq. 12.
As an alternative to its peer ode45 function, ode23 possesses a degree less accuracy while allowing the more efficient processing of stiff ordinary differential equations while still maintaining sufficient accuracy. Using oscillatory shear flow with γ(t) γ 0 sin(ωt) ( Figure 1A), and shear rate as _ γ (t) γ 0 ω cos(ωt) ( Figure 1B), the model predictive capability is assessed. This is accomplished by predicting oscillatory shear flow with both small and large amplitude oscillatory shear flow, using the best-fit parameters. This is shown with SAOS (amplitude sweep), and four sets of LAOS at frequency 1 (rad/ s), at strain amplitudes of 1.0, 5.0, 10.0, and 50.0 (-). Cost functions are computed for each set of LAOS as: the small amplitude oscillatory shear is predicted, and F cost, SAOS is again reduced in accordance with the following expression F cost, SAOS is then calculated for the amplitude sweep, with fit parameters held constant to product a full alternance period.
The expression A = bx is then applied, where A represents and array of two columns: γ 0 sin(ωt) and γ 0 cos(ωt). Then, b is a stress prediction for a period. The first harmonic module, G ′ 1 and G ″ 1 , are then calculated and subsequently compared to each combination of strain amplitude and frequency of SAOS and amplitude and frequency sweeps. This leads to the calculation of the third harmonic moduli, G ′ 3 and G ″ 3 . The various cost functions for the step-ups and step-down fits, SAOS and LAOS predictions are then collectively used to gauge model efficiency.
To fit the tensorial model variants, the following is done; first, initial guesses are provided from the best fit values of the non-tensorial model variants. Next, the steady state and transient parameters are simultaneously fit to a nonnormalized steady state cost function from a collection of stress measurements.

Sequence of Physical Processes
We introduce Sequence of Physical Processes (SPP) with a more in depth discussion of oscillatory shear flow, where the strain is give as γ(t) γ 0 sin(ωt) with γ 0 (-), the strain amplitude, and ω (rad/s) as the frequency of oscillation. The first temporal derivative of strain is the shear rate, given as _ γ (t) γ 0 ω cos (ωt). Figures 1A,B below is the strain vs time, and shear rate vs time; noting that a full period of oscillation at a frequency of 1 (rad/s) is 6.28s.
We must now create the SPP-A matrix with the data directly from the rheometer: the strain; quotient of shear rate over frequency; and stress (all) over a period of oscillatory shear flow at alternance is used to construct the A matrix as follows: . The A matrix is then operated on by calculating the tangent vector, normal vector and binormal vector at each discrete data point from the rheometer over a period of alternance [50][51][52]. Note below that as the strain amplitude gets larger, the shape of the viscous projections in Figure 2B become closer and closer to a perfect diagonal, signifying that the material is getting more and more "liquid-like," or purely viscous, while the shape of the elastic projections in Figure 2A is becoming closer and closer to a circle signifying (again) a more and more "liquid-like" or viscous signal. Figure 2 shows consistency with respect to liquid-like properties.
The tangent: T A _ A , normal: N T _ T , and binormal: B T × N; where the single dot is the dot product, and the x notation is the cross product [50][51][52][53]. Finally, the components of the binormal are used to compute the transient elastic and viscous moduli, With the understanding is that the transient elastic modulus, or the instantaneous 'solid-like' property of the material is equivalent the partial derivative of stress with respect to strain, G ′ t (t) zσ(t) zγ(t) ; and the transient viscous modulus, or the instantaneous "liquid-like" property of the material is equivalent to the partial derivative of stress rate with respect to shear rate, . The Frenet-Serret apparatus applied to the TNB (tangentnormal-binormal) frame is consolidated, and shown below in Table 5, 6. We incorporate the elastic moduli, G ′ t , G ″ t (Pa) and the transient phase angle δ t tan −1 ( radians) to interrogate the material during LAOS flow, correlating these mechanical property metrics with our rheological model metrics of rouleaux (microstructure) [50][51][52]. Moving forward we also make the following substitutions for viscous transient modulus:

Parametric Correlations
Our parametric correlations demonstrated here involve the matrix of all ten donors best t-ESSTV parameter fits, along with the hematocrit (HCT), fibrinogen, (fib.), total cholesterol, triglycerides (tri.), high-density lipoprotein (HDL), and lowdensity lipoprotein (LDL). The Matlab corrcoef command is invoked, whereby the correlation coefficients, and p-values for the matrix of rheological and physiological parameters are calculated. In addition, we show the average, standard deviation, and 85%, 90 % , and 95% confidence intervals are calculated as follows: CI x + /−Z s n √ , where x is the average, s is the standard deviation, n is the total number of donors (10), and Z = 1.44, 1.64, 1.96 for 85%, 90% and 95% respectively. A p-value of less than 0.05 is considered statistically significant [47].

Donor 1 Model Fits and Predictions
Figures 3A-E shows the best fit parameter values of Donor1, using the t-ESSTV model. Figure 1A is the steady state and with steady state microstructure prediction, while Figures 3B,C shows the best fit step-ups from an initial shear rate of shear rate of   The plastic strain rate at the shortest and longest times is equal to the applied rates. While the plastic rate changes monotonically during the step down, it displays both increases and decreases during the step up in applied rate. As expected, with a step-down in shear the microstructure is rebuilding. Recall that the parameters are first fit to the nontensorial version of the model, to obtain a good initial guess. The initial guess is then adjusted by fitting to the tensorial version of the model simultaneously to the steady state and set of three step-up and three step-downs in shear rate experiments. This procedure yielded the fits shown below in Figure 3.
With the determination of the full set of parameters, shown in Table 6, for Donors 1-12 [46], the model is then used to predict the first normal force, N1 shown in Figure 4A, and the effective times of microstructure aggregation and breakdown comparison in Figure 4B. As is shown in Figure 4B the effective time of shear induced structure breakdown is much shorter than the effective time of shear induced aggregation. The physics of rouleaux evolution are the most interesting result of this research. We can now, using components Eq. 4 explain the dynamics of λ shown in Figure 3A. As a result of τ eff,break dominating τ eff,aggr the microstructure, or rouleaux, are driven close to zero for high shear rate. In other words, the effect of shear aggregation, and Brownian motion is not enough to overcome the tendency for the rouleaux to dis-aggregate, with the effect more pronounced as the shear rate is increased. Of course, this is compatible with biological systems whereby it is desirable not to have robust rouleaux formation in all biologic circumstances. This effect is magnified as the shear rates increase moving from left to right across the x-axis, thusly explaining the drop in steady state microstructure values. This explains why the microstructure asymptotically approaches zero at the higher values of shear rate.
Using the t-ESSTV model, a sequence of LAOS and UDLAOS is predicted below, again to demonstrate the efficacy of the t-ESSTV model, shown in Figures 5A,B. With respect to the UDLAOS the strain is now given as γ(t) γ 0 sin(ωt) + γ 0 ωt, and the shear rate is the first-time derivative, _ γ(t) γ 0 ω cos(ωt) + γ 0 ω. Note that the UDLAOS does not change direction yet still undergoes a speed-up and slow-down phase of each period. One could also recast UDLAOS as a linear superposition of constant shear rate and an oscillatory shear component, that when added together is always positive. This is more consistent with the pulsatile flow dynamics of human blood. Note Figure 5B below plots the UDLAOS elastic projection with the oscillatory component of strain only, γ(t) γ 0 sin(ωt), and viscous projection with the oscillatory component of shear only _ γ(t) γ 0 ω cos(ωt). Figure 5A presents the LAOS predictions as a sequence of elastic (left side) and viscous (right side) projections over γ 0 1.0, 5.0, 10, & 50 (−); ω 1.0 (rad/s). The t-ESSTV can qualitatively predict all LAOS, and qualitatively performs better the higher the strain amplitude for all frequencies shown. The elastic projections will approach a straight diagonal line for purely elastic/solid-like conditions, while the viscous projections will approach a straight diagonal line for purely viscous/liquid-like conditions. These two extremes evolve from the following framework: σ(t) Gγ(t) for solids (Hookean approach), relating total stress to strain, where the elastic modulus G is the proportionality constant; and σ(t) η _ γ(t) for liquids (Newtonian viscosity approach), relating total stress to shear rate, where viscosity η is the proportionality constant. The viscous projections show that as the strain amplitude increases at all frequencies the signal become closer to a diagonal line, indicating that the rouleaux are not contributing, because they have broken down, and the blood is much more liquid-like. The shape at the smaller strain amplitudes (at all frequencies) indicates that at these combination of strain amplitudes and frequencies, the rouleaux are relevant and contributing to the total stress signal.  [47,[50][51][52][53].

Stress reconstruction
−′liquid − like′ metric −measure of viscous resistance −instaneous time derivated of stress WRT Shear rate Frontiers in Physics | www.frontiersin.org June 2022 | Volume 10 | Article 889065 Figure 5B presents the UDLAOS elastic and viscous projections together. Where the same general trends are observed. Keeping in mind that with UDLAOS the flow is only in one direction, therefore there is not as much structural breakdown per cycle, as the oscillatory nature is an addition process contributing to structure breakdown. The UDLAOS merely speeds up and slows down, then speeds up, without changing the direction of flow, more consistent with the dynamics of human blood flow. Figures 6A,B below demonstrates the structural evolution differences, as predicted by the t-ESSTV model, using best fit parameter values of Donor 5, over a single period of LAOS and UDLAOS respectively at γ 0 0.5, 1.0, 5.0, 10, 50 & 100 (−) and ω 1.0 (rad/s).
The dynamics of structure breakdown and build up are different. Whereby in general Figure 6A shows that the oscillatory nature of changing direction twice per period causes smaller average values of microstructure averaged over the entire period, while the UDLAOS allows for larger peak microstructure values over each period. This is since LAOS changes direction, and UDLAOS does not change direction of flow during the period. This in turn allows for more time for structure buildup and breakdown to occur during the period, thusly leading to larger microstructure peaks in Figure 6B. The dependence of the microstructure parameter, λ, on the elastic strain is shown in Figure 6C,D. Regions of clear linear dependence are seen at small values of λ, while the dependence is more complex at larger values of λ.
By consolidating all the best fit parameter values from Donor 1-12, shown in Table 5, the parameter averages and standard deviations are computed. Additionally, Table 7 shows the 80%, 85%, and 90% confidence intervals. This allows for a psuedodetermination (albeit with a small number of Donors) of defining left and right limits of mechanically "normal" human blood, over the steady state range of shear rates, shown in Figure 7. Outside of the 90% confidence interval one could medically determine that there is an underlying medical condition requiring additional testing, or analysis. Table 7 below now combines the rheological and physiological parameters of all twelve Donors, along with average, standard deviation, and 80%, 85% and 90% confidence intervals.

Model and Physiological Parameter Correlation
Using the first twelve rows of Table 7 above, the correlation matrix is computed with corrcoef command in Matlab. The results are shown below with the correlation matrix and p-values in Table 8. Correlations and p-values highlighted in bold, red numbers are statistically significant (or close) whereby significance is defined as a p-value less than 0.05. This calculation is also a function of the size of the data set. Italicized numbers in Table 8   Increasing the sample size leads to more accurate and robust correlations. Table 8 shows that hematocrit possesses a strong correlation with zero shear viscosity, infinite shear viscosity and rouleaux viscosity. Hematocrit is significantly correlated with all three viscosity parameters. This makes sense as the hematocrit is proportional to the volume of blood that is taken up by red blood cells. As the volume of red blood cells increases, all three of the viscosity parameters increase. Total cholesterol and highdensity lipoproteins are also correlated with Cross time constant, and elastic modulus of the plasma and individual red blood cells. This is further reflected in a less pronounced, though present correlation between hematocrit and rouleaux elastic modulus. Fibrinogen is negatively correlated (to a small extent) with the Cross time constant and the elastic modulus of the plasma and individual red blood cells. We note that with more Donors, we could strengthen these correlations, and elucidate the ties between rheological and physiological parameters.

Donor 5 Mechanical Properties During LAOS
Below, Supplementary Figure S1 shows the transient moduli, and transient phase angle, G ′ t , G ″ t , and δ t , as well as the t-ESSTV model prediction for λ . The transient moduli are computed directly with the LAOS data, while the prediction of microstructure, λ is calculated with the best fit model parameters for Donor 5 using the t-ESSTV model. Supplementary Figure S2A-C are the contour maps constructed with six sets of LAOS data over γ 0 0.5, 1.0, 5.0, 10, 50, and 100 (−) taken at ω 100 (rad/s). Supplementary Figure S2D is constructed with the analogous model predictions of microstructure for the same strain amplitudes and frequency. The figure is an amalgam of six oscillatory cycles at different strain amplitudes. The different strain amplitude regions are demarcated by letters A-F; while each oscillatory cycle is broken into one-eighth radian segments demarcated by numerals 1-8 in Supplementary Figure S2A,B. These numbers correspond with the cycle number labels shown above in Figures 1, 2, 8.
Supplementary Figure S2A demonstrates a slow evolution of the transient elastic modulus to softer and softer values as the strain amplitude increases from left to right on the x-axis. While the elastic modulus starts at a relative "soft" value but increases in hardness from point 1 to 2, and again from point 5 to 6 Recall that the transient phase angle, δ t tan −1 ( , and takes on values of <0.75 rad for more "solid-like" states, and above this is considered more "liquid-like." This evolution to softer blood comes because of less and less microstructure at the higher strain amplitudes shown in Supplementary Figure S2D. We note here SPP metric δ t lines up with the model prediction of microstructure, λ. This effect is more correlated for the smaller values of strain amplitudes shown in Supplementary Figure S1C,D, and the effect dissipates as the microstructure becomes less and less relevant at higher strain amplitude. This is fascinating because we have shown two different metrics of solid-like behavior simultaneously predicting areas of relative maxima for the LAOS with low strain amplitude indicating that our model can predict areas of high microstructure well, and that this microstructure correlates with more "solid-like" mechanical behavior.
Supplementary Figure S1B shows the transient viscous, or loss modulus over the same range of strain amplitudes. The figure shows a general trend to thinner blood at the higher strain amplitudes, due to the decreasing presence of microstructure or rouleaux. While at the smaller strain amplitudes the rouleaux are still present, at some level, and contributing to an increase in viscosity. This trend is clearly displayed in Supplementary Figure  S1B. More microstructure, λ leads to thickening of the blood, and this is cyclic based on location in the oscillatory shear cycle. For example, from point 2 to 3, and from point 7 to 8 the blood is thickening. From Figure 1B, we see that from point 2 to 3, and again from point 7 to 8 the shear rate is passing through zero, and during this portion of the oscillatory shear cycle the structure is rebuilding and contributes to a thicker blood as there are more rouleaux present, with more resistance to flow.  From this we postulate that at the smallest value of strain amplitude we expect there to be relatively more structure, and in turn this structure will overall contribute to larger values of G ′ t , and G ″ t , aka. overall harder and thicker blood during this LAOS cycle. The numerals represent one-eighth radian segments of a LAOS cycle, and are analogous the numerals in Figures 1, 2, 8, as well as Supplementary Figure S2  In Supplementary Figure S3 we show how the stress depends on the elastic strain and the plastic strain rate, as well as the connections between the rouleaux/microstructure λ and the SPP modulus and viscosity for strain amplitudes of 1, 10, and 100. Even though the strain amplitude increases by a factor of 100 between these tests, the elastic strain decreases with increasing applied amplitude. This reflects the effect of higher shear rates to disrupt the rouleaux/microstructure λ. As a result, the rheology becomes dominated by plastic processes, which leads to the situation where the stress is a single-valued function of the plastic strain rate at the largest amplitude, shown in Supplementary Figure S3J. Similarly, as the amplitude is  increased and the shearing breaks down the rouleaux/ microstructure, λ decreases but also can be described by only the SPP viscosity η ′ t . These results provide support for performing future blood rheology studies using recovery rheology metrics that experimentally distinguish between elastic and plastic strains.
We note here that we can now compare the effect of high hematocrit on the solid-like and liquid-like properties of Donor5 by comparing to the identical figures with the aggregated solid and liquid-like properties of human blood shown in Supplementary Figure S4, below.

Average Human Mechanical Properties During LAOS
The human blood average contour maps shown below were constructed in a similar fashion to Supplementary Figure S1. Supplementary Figure S4 is constructed with the average values of human blood's transient moduli using 12 Donors, borrowing from previously submitted work for comparison [47].
Supplementary Figure S4C was created using the average parameter values of the ten Donors from Table 6. With the average contour maps available one can compare Donor5's mechanical properties in terms of hardness/softness and thickness/thinness with that of the "average". Immediately obvious is the fact that as compared to the average contours below Donor5 has harder, and thicker blood because of the higher-than-average value of hematocrit value of Donor5. We can prognosticate that this could be because that Donor5 has higher hematocrit, and total cholesterol than average here. From previous work it has been shown that the physiological parameters most impacting mechanical properties are hematocrit and total cholesterol [47]. Supplementary Figure  S4C,D show a correlation with the SPP prediction of solid-like behavior and the t-ESSTV model prediction of λ, or rouleaux. The t-ESSTV model has allowed for predictions of microstructure that correlate with predictions of increased "solid-like" behavior as computed by SPP. From this we postulate that at the smallest value of strain amplitude we expect there to be relatively more structure, and in turn this structure will overall contribute to larger values of G ′ t , and G ″ t , aka. overall harder and thicker blood during this LAOS cycle. As the strain amplitude increases to 100 (-) we see the general trend of: less microstructure in Supplementary Figure S5G Figure S2F,K) blood over the LAOS cycle. At the higher strain amplitudes, we can see that microstructure is broken further and further down. We also note that at the higher strain amplitudes in Supplementary Figure S1, S4 there is more red coloring, representing more liquid-like behavior (color-mapping shown in Supplementary Figure S1C, S4C).
Lastly, using the aggregated small amplitude oscillatory shear flow data in terms of G′ and G″, as measured by the ARES-G2, we see that over this range, blood is always more liquid-like than solid-like. Additionally, we superimpose our aggregated model predictions, using the aggregated model parameters, shown in Supplementary Figure S6 below.

DISCUSSION
The full tensor version of the model demonstrated here, t-ESSTV is comprised of thirteen parameters fit simultaneously to steady state with six steps up and down in shear rate. Per previously established methods, three model parameters are fixed with d set to1/2, m set to 3/2, and γ 0,R set to one [7-9, 36, 37, 42, 43, 46]. Figure 3A demonstrates the initial steady state while Figure 3B,C show the three step up and three step down fits, respectively. Meanwhile, Figure 3C,E. reveal the evolution of the structure parameter with variance in shear rate. Table 6 aggregates the best fit model parameters from both steady state and transient data and F cost,ss for steady state, step up and down tests for Donors 1-12. The tensorial nature of the model precludes the derivation of a simple algebraic, steady state solution, as one uses the timeintegration of the ODEs to solve the model equations. As such, the transient and steady state parameters must be fit simultaneously to one set of steady state data along with three sets of step-ups and three sets of step-down transient response tests. Figure 4 lends a picture of the evolving microstructure in detailing the evolution of the structural parameter while Supplementary Figure S5A shows the aggregated SAOS amplitude sweeps. Supplementary Figure S5 as expected suggests a dominant liquid metric, loss modulus G′, at experimental shear rates, as expected of blood. At the experimental conditions, the loss modulus intuitively exceeds the solid-like metric, G″, due to the notable liquidity of the fluid. However, the fit for G″ is notably more adherent to data throughout the range of shear rates while G′ prediction proves far more qualitative. This is a place for future improvements to our model's predictive capability.
In conjunction with the viscoelastic LAOS projections seen in Figure 5, the analysis reveals the relative efficacy of the model at predicting the steady state and transient response of blood under evolving force conditions up to γ 0 50. In this, the tensorial model offers markedly increased predictive efficiency over nontensorial peer models [46]. In inspecting the t-ESSTV parameter fit and F cost values from Table 6, it can be seen that the model enjoys strong performance, with maximally reduced cost values. The low average F cost suggests the utility of t-ESSTV in the effort's objective to characterize TEVP substances such as human blood. The model is shown to predict the SAOS, LAOS and UDLAOS data well, and surely stacks favorably with peer models. This is a benefit, so we can proceed to use the microstructure, or rouleaux, predictions for the LAOS to compare with other objective indications of more solid-like behavior from SPP.
The better the model's predictive capability the more accurate the rheological and physiological parameter correlations will be, as well as the more insight that will be provided. It is this comparison that is enabled by our efforts here. Without a robust model, with strong predictive capability of not only the stress, but also the microstructure (rouleaux), the contour maps of mechanical properties could not be compared to a reliable, objective predictor of areas of more structure. It appears from the contour map analysis that the most correlated regions of microstructure, λ, to δ t , and microstructure to G ′ t happens at the lower strain amplitudes, as expected because at the lower strain amplitudes the structure is not as relevant because the oscillatory shearing as well as the magnitude of the stress here is deleterious to the rouleaux formation.
Lastly, we have compared elastic and plastic metric against the measured stresses in steps up and down in shear rate as well as in LAOS. We have shown that the evolution of the elastic strain and plastic strain rate is non-trivial during steps up and down and can even be non-monotonic. During the largest amplitude LAOS tests, we see in Supplementary Figure S3 that the plastic strain rate and the SPP viscosity are well-suited to a complete description of the stress and microstructural responses. These results provide testable predictions of our model and other models too and lend support for future studies to use recovery rheology protocols that distinguish between recoverable and unrecoverable strains. The recovery measures provide more detailed information than the total strain and total strain rate dependences in isolation, and clearly allow for links to be made between the rheology and the microstructure.

CONCLUSION
In highlighting the development of the t-ESSTV model and its application to the analysis of blood rheology, this effort has demonstrated the model's predictive capability and efficiency, along with further confirming previously established expectations regarding blood thixotropy. Specifically, regarding rouleaux microstructure, the experimental methods have revealed a few rheological and physiological correlations that inform the current understanding of blood mechanics. Specifically, the model results, when processed with Sequence of Physical Processes, it was observed that a notable positive correlation existed between λ, the structure parameter, and known metrics of solid-like behavior such as G′ and δ t . As microstructure aggregates at high values of λ, the apparent increase in solid-like behavior stems from the increasingly relevant thixotropic and viscoelastic effects of aggregated red blood cells. Moving forward our intention is threefold: 1) Continue to expand our database of Donor human blood, branching out into pathological examples when possible; 2) Continue to develop our TEVP modeling capability, with an eye toward introducing a conformation tensor and non-equilibrium thermodynamics (NET) framework to model the effects of the rouleaux more robustly; and 3) Continue to compare elastic, "solid-like" metrics of human blood properties from Sequence of Physical Processes to improved TEVP model predictions of microstructure.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: data.Mendeley.com

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by University of Delaware IRB. The patients/ participants provided their written informed consent to participate in this study.