Novel tensorial Thixo-Visco-Plastic framework for rheological characterization of human blood

Characterizing human blood, a complex material with a spectrum of thixo-elasto-visco-plastic properties, through the development of more effective and efficient models has achieved special interest of late. This effort details the development a new approach, the tensorial-enhanced-Thixo-Visco-Plastic model (t-e-TVP), which integrates elements from the proven Bingham and generalized Maxwell systems to create a more robust framework and subsequently cast into a tensorial format. Here, the elastic and viscoelastic stress contributions from the microstructure are superimposed upon the viscoelastic backbone solution for stress offered by the modified TVP frame. The utility of this novel model is tested against the contemporary tensorial-ethixo-mHAWB (t-ethixo-mHAWB) framework, a similar model with a greater number of parameters, using rheological data of human blood collected on an ARESG2 strain-controlled rheometer. The blood samples are parametrically and statistically analyzed, entailing the comparison of the t-e-TVP and t-ethixo-mHAWB models with their capacity to accurately predict small and large amplitude oscillatory shear as well as unidirectional large amplitude oscillatory shear flow in blood.

www.nature.com/scientificreports/ of lambda values, featuring independent thixotropic evolution timescales, parallel to an isotropic kinematic hardening framework 40 . An analogous effort was accomplished in the modification of the SPTT-Isotropic Kinematic Hardening model ("S" for Saramito's novel plasticity term, "PTT" for the Phain-Thien Tanner viscoelastic model) which acted to combine several approaches to effectively representing the material physics of a TEVP materal 17,41 . These models possessed a tensorial form and possessed 11 to 15 distinct parameters. Concurrent to the developments, the Modified Delaware Thixotropic Model (MDTM) was enhanced with a viscoelastic timescale of the stress response contribution from the component rouleaux, the novel model being dubbed viscoelastic enhanced MDTM (VE-MDTM) 23,31,35,42,43 . This alteration further demonstrated the importance of including the stress response from changes in microstructure. Recent work has also extended toward the tensorial transformation ETV and ESSTV model into tensorial analogues: t-ETV and t-ESSTV respectively 45 .
Despite the relative effectiveness of current generation models in fitting steady state and, most notably, transient rheological data for blood, their relative efficiency in analyzing the material nature of blood does leave something to be desired 5,10 . This effort utilizes elements from the previously established steady-state variants of the Dullaert & Mewis and MDTM models. These are then recast into a dynamic Maxwellian format, allowing for a tensorial representation 17,41 . The development of this new system, dubbed the tensorial-enhanced-Thixo-Visco-Plastic (t-e-TVP) model, is fully described in "Model development" section, where the new method is shown to integrate theories of plasticity to better express the elastic and viscoelastic contributions of the blood rouleaux towards total stress and integrate the full stress tensor 17,40,41 . "Methods" section follows with a description of the experimental protocol relevant to the collection of the experimental samples and a walkthrough of the parametric optimization performed to fit the experimental model to the given data. "Results and discussion" section includes a summary of the experimental results, analyzing the capabilities of the model to predict large amplitude oscillatory shear and uni-directional large amplitude oscillatory shear flow in the circulatory system. The accuracy of the novel t-e-EVP framework in predicting SAOS, LAOS, and UD-LAOS is then compared to that the t-ethixo-mHAWB variant, representing one of the most modern contemporary enhanced rheological modes. The conclusions of the analysis are then enumerated in "Conclusions" section.

Model development
Being of Maxwellian origin, the steady state structure parameter of the tensorial-enhanced-Thixo-Visco-Plastic framework is defined as such: where t r1 and t r2 represent the ratios of rouleaux shear breakdown to Brownian buildup, and the ratio of shear aggregation to Brownian build-up respectively and γ is the fluid shear rate.
(For the tensorial form, γ is given as the second invariant of the rate of strain tensor: γ (1) =γ ij = ∇ν + (∇ν) T ). For this particular application of the steady-state structure parameter equation, the power law of shear aggregation, usually defined as d, is set to ½ per the findings of previously published literature on the subject of the blood medium 8,37,45 . As mentioned above, the steady-state thixotropy value is constrained between 0 and 1, with the former indicating the presence of no bonds between a given RBC and its neighbor and the latter signifying maximum structure agglomeration 14,15,21 .
To capture the elastic and viscoelastic stress contributions from the blood's component rouleaux, we adhere to the approach of Horner et al. in integrating the following equation into the model, where the structural evolution manifests as a function of blood's shear rate 3,8,37,45,46 : where τ is the total rouleaux agglomeration constant and, as before, d is set to ½ in deference to blood's unique properties 8,37,45 . Equation (2) contains three separate components relevant to the thixotropic breakdown and formation of Rouleaux: 1) Shear breakdown proportional to shear rate and the amount of remaining structure, (2) Shear buildup characterized by (1 − ) ; and (3) Brownian aggregation also defined by (1 − ) 3, 46 . In capturing all three of these rouleaux phenomena, the model's capability for concrete analysis is enhanced. Further, the generalized Maxwellian framework is evolved for use in defining the steady-state yield stress of the rouleaux, producing: where σ y,0 represents the yield stress, m represents an additional power law parameter, set to 3/2 in this work, while η ST and η ∞ represent the structural viscosity and infinite viscosity respectively. The t-e-TVP's fundamental description of the stress contributed by individual RBC's possible deformations manifests as a modification and combination of the Kelvin-Voigt: and generalized Maxwellian: yielding: www.nature.com/scientificreports/ where G represents the elastic modulus. To better enhance the stability of the model equations, a tensorial approach is adopted in describing the individual xx, yy, zz, and yx axis components, the full Eqns. (7 -10) shown below 3,8 : With the final introduction of the tensorial components, the description of the novel model is complete. As can be seen the model is fully capable of predicting the first normal stress difference, N1 as σ xx −σ yy . Table 1 contains a prompt summarization of the tensorial-enhanced-Thixo-Visco-Plastic model discussed here.
Being counted among one of the most useful contemporary models the tensorial ethixo-mHAWB, an enhanced tensorial variant of the original mHAWB model, provides a basis of comparison against for the framework developed above 3,44 . The t-ethixo-mHAWB approach requires 13 parameters, 10 of which must be fitted via the procedure discussed in "Methods" section. The model's representation of viscoelastic stress contributions is shown below in Tables 2 and 3.
In utilizing the t-ethixo-mHAWB model for comparison, a clearer understanding of the relative efficacy of the new t-e-TVP model can be developed. The upcoming "Methods" section will detail the experimental procedure relevant to the procurement of the two analyzed blood samples and, subsequently, the methods by which the collected data was fitted to the two subject models.

Methods
Experimental protocol. The prosecution of this effort's blood collection and rheological characterization protocols accord with existing precedent for blood rheology experiments, as previously enumerated by Horner et al 3,5,8,46 . The blood draw procedure was compliant with the expectations of the University of Delaware's institutional review board, with the blood being drawn from the seated patients' antecubital veins 3,5,8 . "The human blood samples were collected by a licensed practitioner at the Nurse Managed Primary Care Center located at University of Delaware STAR campus in compliance with, and approved by UD's Institutional Review Board (Study Number 767478-2) 3-5 . " The methods were performed in accordance with our IRB's relevant guidelines and regulations. Additionally, all human participants involved gave informed consent via face-to-face discussion, and signed documentation verifying such. At the time of the draw, the subjects had been confirmed to be free of any deleterious health conditions and had undergone fasting 8-10 h prior to the sample retrieval. 6 ml of retrieved blood was mixed with 1.8 mg/ml ethylenediaminetetraacetic acid within a vacutainer tube while 9 ml of blood was set aside for a comprehensive blood count, lipid panel, and fibrinogen activity test. The results of the latter testing can be found in Table 4 below 3,5,8 . A TA Instruments ARES-G2 strain-controlled rheometer with a double wall couette was utilized for all experimental measurements 3,5,8 . The specific dimensions and measurable range of the couette geometry are identical to those of previous work 44,46,47 . No later than 60 min from collection the blood samples were introduced to the rheometer for testing, with all analysis occurring within 4 h from withdrawal. Through the tests the temperature was maintained at 37.0 °C via a Peltier temperature controller. Shear rate was consistently maintained below 1000 s −1 and each test was initiated with an applied preshear of 300 s −1 for 30 s to eliminate lingering memory effects on microstructure 3,5,8 .
Experimental results are shown for steady state, step up/down in shear rate tests, amplitude and frequency sweeps, large amplitude oscillatory shear (LAOS) and unidirectional large amplitude oscillatory shear flow  The representative function for both the amplitude and frequency sweep is as such:  www.nature.com/scientificreports/ where the function inputs are strain amplitude γ 0 , time t, and frequency ω . For the amplitude and frequency sweeps, the running index n is customarily set to 1, though the upper summation limit is infinite. Additionally, the primary metric of success for model predictions is found in the ratio of the third to first harmonic: with the subscript on the elastic moduli representing the harmonic order 3,8,37,39 . Used to quantify the efficiency of a given model, Eq. (16) juxtaposes the relative intensity of the third harmonic to the first harmonic. This analysis of model effectiveness is shown on the models' amplitude and frequency sweep predictions in Figs. 3 and 4 3 where y i and f i represent the steady state stress data and model prediction respectively. The transient parameters are then experimentally fit to four tiers of shear rate step up and four tiers of shear rate step down with parallel tempering, acting to minimize the cost function: where M and N are the total number of transient step up/down tiers in shear rate testing and number of datapoints respectively. The particulars of the experimental fitting procedure are like those of Armstrong at Tussing 35 .
The rheological dataset predictions can be seen in Fig. 2 below, with LAOS and UD-LAOS cost functions calculated as such: Amplitude and frequency sweep analysis is conducted with the cost function described as: where n represents the number of datapoint. F cost,Sweep is produced for both amplitude and frequency sweeps, while the model parameters are held constant to predict the full alternating period. This utilizes an array with A = bx where x consists of two elements: γ 0 sin(ωt) and γ 0 cos(ωt) while b is the period's stress prediction 37,39 . The third harmonic moduli, G ′ 3 and G ′′ 3 , are also derived as shown in Eq. (20). All of the parameter optimization for both models was conducted with Matlab, version R2021a, with a stochastic, global, parallel tempering algorithm. The data for the figures was obtained with Matlab, and the final figures were constructed with Origin Graphing and Analysis software version 2021b 31 .
The analysis of the LAOS and UD-LAOS integrate the experimentally derived F cost,LAOS over three data period, with the final period predictions and data demonstrated in the Lissajous-Bowditch elastic and viscous projections. The analysis procedure utilizes both the Akaike Information Criteria (AIC) and the Bayes Information Criteria (BIC) (with steady-state and transient step optimization) to juxtapose additional parameters against model accuracy, penalizing the addition of extra parameters 32,35,43,48 . The AIC and BIC residual sum of squares (RSS) is calculated as such: where y i and f i are the data and model prediction respectively. For AIC and BIC, the steady state cost function and RSS is normalized by 32,35,43,48 :

Results and discussion
t-e-TVP model. The t-e-Thixo-Visco-Plastic framework possesses 9 total parameters, with two being fixed with reference values, five being fit with steady-state analysis, and the final two resulting from the transient step up/down shear rate analysis. The two fixed values, d, and m are set to ½ and 3/2 respectively, as previous literature recommends for analysis of blood 3,8,37,43,46 . Due to the tensorial nature of the novel model equations, a trivial, algebraic steady-state solution is unnecessary. As such, the transient and steady-state parameters are simultaneously fit to one set of steady-state data and six sets of transient step down/up shear rate data. Figure 1a details the steady state fit while Fig. 1b,d show the three iterations each of the step down/up transient shear rate fits. From these two plots, it can be deduced that the model is generally sufficient in predicting the stress under evolving shear rate of the experimentation. In addition, Fig. 1c,e showcase the corresponding structure parameter's evolution as it is subjected to steps down and up in applied shear rate. The optimized fit model parameters from the steady-state and transient analysis, as well as the F cost values for step down/up and amplitude/frequency sweeps, can be seen in Table 5. Table 5 Figure 2c,e showcase the structure parameter's corresponding evolution as it steps down and up in shear rate. The t-ethixo mHAWB best fit model parameters from the steady-state and transient analysis, AIC, BIC, RSS, and the F cost values for step down/up, amplitude sweeps and frequency sweeps, can also be seen in Table 5. F cost, is a dimensionless value due to having been normalized at each data point. Lastly, Fig. 4a Table 6 will reveal that there exists a tangible Table 5. Best fit model parameters, and F cost of fits for Donor 1 with t-e-TVP and t-ethixo-mHAWB models. Highlight legend: yellow-steady state fit; green-fixed; orange (salmon column) t-e-TVP step up/down fits; orange (blue column)-t-ethixo-mHAWB step up/down fits. (Donor 1, Dataset 1) 49 .

DONOR1
Par. interpretation of the elastic and viscous projection of LAOS and UDLAOS demonstrate the predominance of viscous effects, as is expected from the blood medium 8 . This interpretation is derived from the legends contained in Fig. S5 of the attached supplemental materials. Nonetheless, in this analysis it is important to remain cognizant that the demonstrated data is only that of Donor1, with figures for Donor 2 located in the supplemental materials. Figures 3a,b and 4a,b shows that both models perform nearly the same in prediction the storage and loss modulus, for the amplitude and frequency sweeps respectively. We note here that blood over this range of strain amplitudes and frequency always shows a larger value of loss modulus (liquid-like metric), G" than (solidlike metric) G' , as is expected, as human blood over this experimental range in shear rate, and temperature is more a liquid, than a solid. For both Figs. 3a,b and 4a,b the G" values are predicted more accurately, than is G' for both models, however the G' trends are qualitatively predicted by both models. This is because human blood is challenging to work with since at most of the combinations of strain amplitude and frequencies the blood is much more 'liquid-like' than 'solid-like' . This is reflected in the actual data and the model predictions. Furthermore, with respect to the amplitude sweep, to the left of values of strain amplitude: γ 0 = 1(−), there is clearly some vestigial structure (rouleaux) present, while after this value there is not. On top of that, at the frequency of the amplitude sweep, ω = 12.566 (rad/s), and γ 0 ≤ 1(−), values less than this is in the vicinity of the linear region of SAOS. The 'doublepeak' shown of the I 3 /I 1 trend, most likely have to do with the fact that in this region there are rouleaux evolving:  Table 6 quantitatively demonstrates the same. In comparison of Figs. 3c and 4c (UDLAOS), and Figs. 3d and 4d (LAOS) we qualitatively see that both models' predictive capability of UDLAOS and LAOS is almost identical until the point of strain amplitude equal or exceeding 10(−). At values of γ 0 ≥ 10(−) the t-e-TVP is no longer able to quantitatively predict accurate stress values, while the t-ethixo-mHAWB can. This result is seen quantitatively in Table 6 and corroborated with Donor2 fits and predictions shown in the Supplemental Material. (Analogous fits and predictions for Donor 2 using both models are shown in the Supplemental Material).

Conclusions
This effort has manifested in the development of a novel approach to the problems posed by TEVP fluids such as blood in the creation of a model derived from the modification and subsequent modification of the Kelvin-Voigt and generalized Maxwell frameworks. In applying precedent methods established in the work of Armstrong et al., Varchanis et al., and Wei et al., both modified models were enhanced with a full tensor, granting them relatively superior predictive performance under LAOS and UDLAOS conditions for both Donors 1 and 2 (see supplemental materials below) [39][40][41] . Moreover, this improvement merited the collection of no additional parameters. As has been found in previous literature, the improvements offered by the full tensor are marked 45 . While the predictive capabilities of the t-e-TVP model are not precisely equivalent to that of t-ethixo-mHAWB, analysis   www.nature.com/scientificreports/ of both Datasets 1 and 2 indicated that the t-e-TVP model offers the benefit of drastically reducing the number of fitted parameters, possessing 9 total as opposed to its peer's 13 total. This manifests as t-e-TVP's considerably lower AIC and BIC values, offsetting the minor accuracy advantage of t-ethixo-mHAWB. Future efforts will likely pertain to the further optimization of the novel model's performance through further modification as despite its relatively simple framework, necessitating few parameters, it offers considerable predictive potential. Possible future models will retain the full tensorial approach, with any considered parameter modification or additions grounded in the concrete mechanics of viscoelastic phenomena and microstructure evolution.