Non-identifiability of parameters for a class of shear-thinning rheological models, with implications for haematological fluid dynamics

Choosing a suitable model and determining its associated parameters from fitting to experimental data is fundamental for many problems in biomechanics. Models of shear-thinning complex fluids, dating from the work of Bird, Carreau, Cross and Yasuda, have been applied in highly-cited computational studies of heamodynamics for several decades. In this manuscript we revisit these models, first to highlight a degree of uncertainty in the naming conventions in the literature, but more importantly to address the problem of inferring model parameters by fitting rheology experiments. By refitting published data, and also by simulation, we find large, flat regions in likelihood surfaces that yield families of parameter sets which fit the data equally well. Despite having almost indistinguishable fits to experimental data these varying parameter sets can predict very different flow profiles, and as such these parameters cannot be used to draw conclusions about physical properties of the fluids, such as zero-shear viscosity or relaxation time of the fluid, or indeed flow behaviours. We verify that these features are not a consequence of the experimental data sets through simulations; by sampling points from the rheological models and adding a small amount of noise we create a synthetic data set which reveals that the problem of parameter identifiability is intrinsic to these models.


Introduction
Many complex fluids exhibit shear rate-dependent viscosity; suspensions, in particular fluids of biological importance such as blood, and biological polymers, such as mucus, are typically shear-thinning (pseudo-plastic), i.e. their viscosity reduces with increasing shear rate. A number of models have been proposed for this behaviour and studied intensively; we will focus on a class of models which relate the shear viscosity to shear rate via nonlinear algebraic equations, in particular the formulations of Cross, Bird, Carreau and Yasuda, and their subsequent application to blood rheology. We will address two significant issues -first, an inconsistency in the literature regarding naming of models, and more importantly, some significant difficulties which appear in determining model parameters through least squares fitting. Since there are major (and unexpected) differences in parameter identifiability between subtly different models, unambiguous naming will turn out to be very important.
To set the scene we will briefly review the key models.
The earliest model of Ostwald (1925) and de Waele (1923) is based on a power-law dependence of viscosity on shear rate; limitations of this simple model include is its singularity at zero shear rate and inability to capture high shear rate dependency when compared to empirical data. For these reasons we will not consider the model further. These discrepancies were addressed by Cross (1965) who postulated a four parameter, constitutive relationship: where µ is the effective viscosity of the fluid as a function of shear rateγ, the parameters µ 0 and µ ∞ are the zero and infinite limit shear viscosities respectively, λ is a constant with dimensions of time, and n is the power-law index (equation (1) is presented in a slightly different form from Cross (1965) though is functionally identical). The model has finite, non-zero viscosity, at both zero and infinite shear rate limits. By contrast, the three parameter model of Carreau (1972) provides a finite viscosity at zero shear rate, and zero viscosity in the infinite shear rate limit, The original paper uses the parameter S = (1−n)/2. In a study of polystyrene fluids, Yasuda (1979) modified this formulation to include a further parameter a to describe better the low shear to power-law transition region: This model has four free parameters (µ 0 , λ, n, a), and implies a zero viscosity limit as shear rate tends to infinity. Perhaps surprisingly, the canonical text of Bird et al. (1987) (while citing the same sources as above) defines a different model as the 'Carreau-Yasuda' model µ Equation (4) differs from both Carreau's and Yasuda's models through the inclusion of an infinite shear rate viscosity parameter µ ∞ (in the manner of Cross (1965)), amounting to five parameters. Bird et al.'s five-parameter "Carreau-Yasuda" model (4) has been used for blood flow modelling in key papers by Perktold & Rappitsch (1995), Gijsen et al. (1999) (who referred to Bird et al. and also termed it Carreau-Yasuda) and Leuprecht & Perktold (2001) (who referred to it as a modified Cross model). Equation (4) can be viewed as a hybrid of Carreau, Yasuda and Cross' contributions, which perhaps explains the proliferation of terminology. Indeed for suitably-chosen parameters, equation (4) can be reduced to each of the preceding models.
Most rheological literature takes the slightly different approach of fitting the model to data on a log-log scale. In the framework of maximum likelihood estimation, this approach corresponds to the assumption of lognormal error. Mathematically, we may define, Data from Skalak et al. (1981) (extracted from Ballyk et al. (1994) using GRABIT, Doke (2016)) will be used in what follows because of its excellent coverage in shear rate, from 0.1 s −1 to 500 s −1 .

Results
Each model will be considered in turn, starting with the model possessing the fewest free parameters. The units of µ 0 , µ ∞ will be centipoise, λ will be seconds, and n, a are dimensionless.

Carreau-1972 three-parameter model
Performing a maximum likelihood fit of the Carreau-1972 model with normal errors, i.e. by minimizing L, immediately yields difficulties with parameter identification. Constrained optimization (Matlab fmincon) consistently finds a solution at the boundary of the parameter space (either the maximum value of µ 0 or the maximum value of λ depending on the limits chosen, however with n = 0.483 consistently. The reason for this behaviour is evident in Figure 1a; the likelihood surface L(µ 0 , λ, * ) exhibits an extended 'ridge'. Taking an upper bound of µ 0 < 350 yields the parameter estimate µ 0 = 301, λ = 200, n = 0.483, depicted with a blue star; the cost function value is L * = 25.8. To show how indeterminate this fit is, we will examine an arbitrarily chosen 'alternative' parameter tuple of µ 0 = 211, λ = 100 and fitted optimum n + = 0.482 (red plus), which has a very similar cost function value of L + = 25.9. For any practical purpose, the fits are identical, as shown in Figure 1b,c. The data do not therefore reliably constrain the parameters µ 0 or λ. It is also of note that either parameter set fits the data well with these parameters up to approximatelyγ = 10 s −1 , however they both perform rather badly for higher values of shear rate.
One may ask whether the more traditional approach of fitting to the loglog plot, i.e. by taking lognormal error and minimizing , might work better. The results of this process are shown in Figure 2. The higher shear rate region (10-500 s −1 ) is fitted much better, however the indeterminacy issue is still present. The best fit tuple found by fmincon (with the same bounds) is µ ♦ 0 = 137, λ ♦ = 200, n ♦ = 0.635 which has cost function ♦ = 0.732. A manually and arbitrarily-chosen tuple µ × 0 = 107, λ × = 100, n × = 0.634, plotted as a red cross, yields a very similar cost function value of × = 0.733. Again the flow curves corresponding to each parameter set are essentially identical (Figure 2a,b). While the fit is arguably better than for Figure 1, the parameters µ 0 , λ are again indeterminate from the data. The same issue occurs for several other experimental blood rheology data sets, see A).
Is this problem a consequence of the data available, or is it intrinsic to the Carreau-1972 model? We generated a synthetic data set by choosing parameter values and simulating an experimental series of 50 samples, taken over an extremely wide range of shear rates (10 −3 -10 3 s −1 ) and with a very small addition of lognormal noise (with standard deviation 0.02). The synthetic data are shown in Figure 3a, with fitting results in Figures 3b,c. Fitting over the full range in shear rate (Figure 3b) reveals that a local minimum is now evident, although with again a rather elongated basin. However restricting the range in shear rate at the lower end to 10 −2 -10 3 s −1 (Figure 3c), the basin is extended in a similar way to the real data fit of Figure 2. While this shows that it is possible to alleviate somewhat the indeterminacy, for blood rheology it requires a range of low shear rate data that may not be possible to achieve in practice. We are not aware of measurements of this type being available in the literature. Cost function, optimized over n for each (µ 0 , λ) tuple for synthetic data generated over the shear rate range (10 −3 -10 3 s −1 ); (b) Cost function for 50 synthetic data points generated over the shear rate range (10 −2 -10 3 s −1 ).

Cross-1965
As described above, the Cross-1965 model involves the parameters µ 0 , λ, n, and, in addition, an infinite shear rate viscosity µ ∞ . To facilitate comparison with the Carreau-1972 model we will initially set µ ∞ = 0, which we will refer to as the 'Cross-Zero' model, The result of fitting this model with lognormal error to the data of Skalak et al. is shown in Figure 4a,b -while the fit is excellent, the minimum of the cost function again appears at the boundary of the domain. Extending the bounds of the search space has similar effects to the Carreau-1972 model. The reason for this can be seen by inspecting equation (8): for sufficiently large λ and non-zero shear rate, the constitutive law can be approximated by, yielding an infinite family of approximately equivalent parameterizations for which µ 0 /λ 1−n is constant.
Having established that this simplified version of the Cross model is also affected by parameter indeterminancy, we turn our attention to the full fourparameter Cross-1965 model ( Figure 5). While the fit (Figure 5a,b) is rather better than both the Carreau-1972 and Cross-Zero models, particularly for larger shear rates, a similar parameter indeterminacy occurs as for Carreau-1972 and Cross-Zero (Figure 5c).

Yasuda-1979
The Yasuda-1979 model differs from Carreau-1972 only through an additional index parameter. An interesting effect of including this parameter is that a local minimum is now found interior to the search domain (Figure 6), at µ 0 = 106, λ = 98.1, n = 0.635, a = 7.61. Nevertheless, there is still a very elongated ridge in parameter space and associated uncertainty.

Discussion
This paper considered the identification of model parameters from experimental data for various cases of what we have termed the Bird-Cross-Carreau-Yasuda class of steady shear-thinning rheological models, specifically applied to blood data. Given that all of the models considered exhibited significant uncertainty regarding parameter values -and in the case of the Carreau-1972 andCross-1965 models the optimum value depended entirely on the specification of the search domain -it is clear that it is necessary to be cautious regarding the physical interpretation of the parameters derived from such a fit. While the flow index n was very consistent, the parameters µ 0 and λ are indeterminate and therefore cannot be used to draw conclusions about the zero shear viscosity or relaxation time of the fluid. One may ask whether this parameter indeterminacy actually matters for flow simulation. After all, if one has an accurate model of the response of the fluid to a range of shear rates, why would the individual values of the parameters used to produce this curve matter? To provide insight into this question, we computed pressure-driven axisymmetric pipe flow with the Carreau-1972 andBCCY-1987 models with each of the 'best fit' and 'alternative fit' parameter choices. The results are shown in Figure 8. In all cases the pressure gradient was chosen as 10 dyn/cm. For each case there is a significant relative difference between the flow profiles for each parameter fit. Parameter indeterminacy may therefore significantly affect flow predictions, particularly for flows involving low shear rates.
The rheology of shear-thinning fluids, and indeed the specific field of blood rheology, are much wider-ranging than the class of models and steadyshear experiments we have considered here. For recent review, see Anand & Rajagopal (2017) and examples of recent on time-dependent flows and extensional rheometry, see references Apostolidis et al. (2015) and Kolbasov et al. (2016). The importance of the Bird-Cross-Carreau-Yasuda class of models is underscored by the fact that they have formed part of major works such as the highly cited papers of Perktold & Rappitsch (1995) and Gijsen et al. (1999). As such the matter of parameter identification for these models is important to address. Our investigation has shown that parameter fitting (c) (d) Figure 8: Velocity profiles computed for pipe flow due to pressure gradient 10 dyn/cm for the best and alternative parameter fits for (a,b) Carreau-1972, (c,d) BCCY-1984. (a,c) pipe radius 0.1 cm and (b,d) pipe radius 1 cm.
for this class of models is indeterminate for several key blood rheology data sets in the literature. Moreover, this is still an issue even for simulated high accuracy data over a very wide range of parameter sets. The implications of this indeterminacy are that parameter values for µ 0 and λ in particular cannot be physically interpreted, and predictions from pipe flow models may also be subject to uncertainty. In future studies -not just involving the BCCY class of models -it will be important to assess parameter sensitivity in order to have confidence in model predictions.

Data accessibility
All data and code for generating the figures in this report can be accessed in the GitLab repository: https://gitlab.com/meuriggallagher/nonnewtonianparamident