# Possible counter-intuitive impact of local vaccine mandates for vaccine-
preventable infectious diseases

Maddalena Donà and Pieter Trapman
Faculty of Science and Engineering, University of Groningen, NL

###### Abstract

We model the impact of local vaccine mandates on the spread of vaccine-
preventable infectious diseases, which in the absence of vaccines will mainly
affect children. Examples of such diseases are measles, rubella, mumps and
pertussis. To model the spread of the pathogen, we use a stochastic SIR
(Susceptible, Infectious, Recovered) model with two levels of mixing in a
closed population, often referred to as the household model. In this model
individuals make local contacts within a specific small subgroup of the
population (e.g. within a household or a school class), while they also make
global contacts with random people in the population at a much lower rate than
the rate of local contacts.

We consider what happens if schools are given freedom to impose vaccine
mandates on all of their pupils, except for the pupils that are exempt from
vaccination because of medical reasons. We investigate how such a mandate
affects the probability of an outbreak of a disease and the probability that a
pupil that is medically exempt from vaccination, gets infected during an
outbreak. We show that if the population vaccine coverage is close to the
herd-immunity level then both probabilities may increase if local vaccine
mandates are implemented. This is caused by unvaccinated pupils moving to
schools without mandates.

Keywords— Stochastic SIR epidemic, household model, vaccine preventable
infectious diseases, vaccine mandates, branching processes

##  1 Introduction

###  1.1 Background

Measles is a viral disease that yearly kills over a hundred thousand people
worldwide [17]. It is one of the most contagious diseases known, with a basic
reproduction number $R_{0}$ estimated to be between 12 and 18 [15] (see
Section 1.2 and [7] for a definition of $R_{0}$). Because measles is extremely
contagious and immunity after infection is typically life-long, in the absence
of vaccination, most people will be infected before reaching adulthood, hence
measles is often called a childhood disease.

In some countries measles and other vaccine-preventable infectious diseases,
such as rubella, mumps and pertussis, are still endemic. That is, they are
always present and cause regular peaks of prevalence. In many other countries,
measles is largely absent because of successful vaccination campaigns [17,
16], which caused herd-immunity in the population [7, p69]. We focus on those
countries. In several of the countries in which there have not been
considerable outbreaks of measles for decades, the vaccine coverage has been
declining, and there is fear that the “pool of susceptibles” will grow big
enough to loose herd-immunity and allow for new waves of infection.
Governments and public health authorities try to tackle this by campaigns to
promote vaccine uptake and sometimes by legislation which makes some vaccines
compulsory for children, e.g. in Italy [13], or considering giving
schools/daycare centers the freedom to use a vaccine mandate, as was
considered in the Netherlands [8].

In this paper we consider the case in which vaccination is not mandatory for
everybody and the decision of requiring the vaccine is left to the schools as
was proposed in the Netherlands. In this situation, children who cannot get a
vaccine for medical reasons can freely choose where to attend school, whether
it has a mandate or not. We refer to those children as medically exempt
children. So, in this situation unvaccinated children that are not medically
exempt all go to schools without a mandate, and in those schools the fraction
of unvaccinated children may be higher than if vaccine mandates are not
allowed, which (as we show) may lead to higher risk of outbreaks of the
vaccine-preventable diseases.

It is likely that medically exempt children go to schools that require
vaccination, in order to obtain a safer environment. We are interested in the
question whether and under which conditions local vaccine mandates will lead
to increased risk of infection for unvaccinated children that are medically
exempt.

###  1.2 Model

In this paper we use a mathematical model to analyze the spread of measles
among elementary school-aged children, where we take into account that most
spread of the pathogen takes place within classes. Although we take measles as
the pathogen of interest, the analysis applies equally to other highly
infectious vaccine preventable diseases. We use a stochastic SIR (Susceptible,
Infectious, Recovered) epidemic model in which the relevant population is
closed and partitioned in school classes. So, the total number $N$ of children
is fixed and we do not consider migration, births or deaths and we assume that
the time scale of an outbreak is such that we may ignore children moving to
different classes. We define $n_{c}$ as the number of children per class,
which we, for convenience of exposition, consider fixed, and $m_{c}$ as the
number of classes. So $N=n_{c}m_{c}$. In Appendix A we show that the
mathematical analysis is easily expanded to allow for variation in the class
sizes. We remark that strictly speaking we analyse our model for a sequence of
populations, where $N\to\infty$.

In SIR models each individual can be susceptible to the disease, infectious or
recovered and immune. A contact trough which the virus is spread is called an
infectious contact. An infectious individual remains so for a random length of
time that is called the infectious period, during which they are able to
infect other individuals in the population. In our analysis we only focus on
the probability that an introduction of a virus may lead to a large outbreak
and on the probability that a given child gets infected. We do not consider
the dynamics of such an outbreak in real time. Because of this, possible
exposed/latent periods, between a child getting infected and becoming
infectious [1] have no impact on our quantities of interest.

In many epidemic models it is assumed that the population is homogeneous,
which means that the individuals mix uniformly. However, in real life this is
rarely the case. Here we introduce heterogeneity by allowing individuals to
have local (within class) and global contacts (between classes) with different
probabilities. We assume that, conditioned on an individual being infected,
whether or not infectious contacts are made with different individuals are
independent. So, we implicitly assume that the infectious period is not
random. Obviously, the infectious period for measles is random. However, the
variance of time between the start of the infectious period and the onset of
rash (which is close to the period that a child stays socially active and is
likely to infect other children) is very small [10]. This makes that we can
simplify the analysis and assume that contacts between individuals occur
independently [12]. In our model an infectious individual makes an infectious
contact with each classmate with probability $p_{L}$ and with other
individuals in the population with probability $p_{G}$. This model is also
referred to as the household model, see e.g. [3] or [1, Section 6.3], since
people are divided into small groups (such as workplaces, households or
classrooms) where transmission of the virus is more frequent. The global
transmission probability is typically much smaller than the local probability,
and in our case we assume that $p_{G}$ is asymptotically inversely
proportional with the population size. That is, if we consider a sequence of
epidemic models indexed by the population size $N$, then there exists a
constant $\lambda_{G}\in(0,\infty)$ such that
$\lim_{N\rightarrow\infty}Np_{G}=\lambda_{G}$.

The basic reproduction number $R_{0}$ is the expected number of infections
caused by a “typical” infected individual infected during the initial stages
of an epidemic in a population that is mostly susceptible [7]. (So, $R_{0}$
refers to a situation in which there is no vaccination.) In homogeneously
mixing populations defining $R_{0}$ is straightforward. Although there is some
issue with deciding what a typical infected individual is, it is also possible
to define $R_{0}$ and other reproduction numbers for household models [14, 4].
In Section 2 we deduce the model parameter $\lambda_{G}$ from estimates of
$R_{0}$ from literature.

For measles, the basic reproduction number $R_{0}$ is estimated to be between
12 and 18 [15, 10], which makes the virus one of the most contagious viruses
in the world. We simplify our model by using the approximation $p_{L}=1$,
which means that if one of the unvaccinated children in a class gets the
disease, then they will infect all of their unvaccinated classmates. Setting
$p_{L}<1$, would complicate the exposition considerably, while the qualitative
behavior of the model stays similar, as long as $p_{L}$ is large enough (say
larger that $4$ divided by the class size, when the probability that entire
class gets infected is roughly $98\%$), see Appendix B.

Because $p_{L}=1$ a child and their non-vaccinated classmates together can be
seen as a “super individual”, of which the susceptibility and infectivity is
proportional to the number of unvaccinated children in the class. We assume
that the latent period of measles is relatively long, so that we can safely
assume that all secondary infections in a class are directly caused by the
initially infected child in the class, and not by another secondary infected
in the population. Lastly, we define a time unit as the fixed infectious
period [1, p4].

To analyze the impact of vaccination (or the lack thereof), we assume that the
vaccine against measles is perfect, i.e., once vaccinated a person can neither
get infected nor infect other people. So, vaccinated people are effectively
omitted from the model. A fraction $v$ of the population is vaccinated; we
call $N_{u}(v)=(1-v)N$ the number of unvaccinated children, which are the
initially susceptibles in our model. Another way to introduce vaccines is to
say that each child is independently vaccinated with probability $v$. The two
assumptions lead to slightly different models, but since they are
asymptotically equivalent, we use one or the other depending on what is more
convenient in the situation. For example, when it is useful to think about
$N_{u}$ as a fixed number, we will use the first assumption; if we want the
number of initially susceptible children in different classes to be
independent, then the second model is better.

In order to protect children at a school, the school may (if allowed by law)
institute a local vaccine mandate, which require all pupils that are not
medically exempt to be vaccinated. In our model we denote the fraction of
classes with such a mandate by $\pi$.

Further recall that we require $\pi\leq v$ for reasons of consistency.

##  2 Analysis

To evaluate the global infection rate $\lambda_{G}$ for given $R_{0}$, we
consider a situation with no vaccines nor mandates, i.e., both $v$ and $\pi$,
the fraction of classes with a mandate, are set to zero. Recall that classes
are all of fixed size $n_{c}$. We distinguish between children that are
infected through a global contact and children that are infected within a
class through a local contact. The first group of children are the primary
cases or type-one individuals, that are the first individuals of their classes
to get infected; each primary case will then infect all the $n_{c}-1$ other
(unvaccinated) children in their class, which we refer to as secondary cases
or type-two individuals. We assume that both types make global contacts (and
cause primary cases) at rate $\lambda_{G}$, but only primary cases cause
secondary cases. For this household model, it is known that the epidemic
threshold parameter $R_{0}$ is the dominant eigenvalue of the next generation
matrix

$A=\begin{bmatrix}\lambda_{G}&n_{c}-1\\\ \lambda_{G}&0\end{bmatrix}$ (1)

where element $A_{i,j}$ is given by the mean number of type-$j$ individuals
that a type-$i$ individual infects during their infectious period ([5, 3]). So
$R_{0}$ solves

$-(\lambda_{G}-R_{0})R_{0}-\lambda_{G}(n_{c}-1)=0.$

So,

$\lambda_{G}=\frac{R_{0}^{2}}{R_{0}+n_{c}-1}$ (2)

If $R_{0}=15$ and $n_{c}=25$, this leads to $\lambda_{G}\approx 5.8.$

Note that in [5], the dimension of the next generation matrix is $n_{H}\times
n_{H}$, where $n_{H}$ is the maximal household size. In our case it is enough
to consider a $2\times 2$ matrix because we only need to divide individuals
into primary and secondary cases. In Appendix B, when dealing with $p_{L}<1$,
we need higher dimensional matrices.

The dominant eigenvalue of the matrix $A$ is sometimes called the individual
reproduction number and is denoted by $R_{I}$ [4]. Assuming $p_{L}=1$, the
reproduction numbers $R_{0}$ and $R_{I}$ are the same because all secondary
cases are infected by primary ones only.

For the analysis of the impact of vaccination in the population, we call the
fraction of the population that is vaccinated $v$, and we define $\pi$ as the
fraction of classes with a vaccine mandate. As mentioned before, there are in
total $m_{c}$ classes, each with $n_{c}$ students. So, $N=m_{c}n_{c}$. The
number of classes without mandate is $m_{c}(1-\pi)$ so the maximum possible
number of unvaccinated children is $m_{c}n_{c}(1-\pi)=N(1-\pi)$. Therefore,
the number of unvaccinated children $N(1-v)$ has to be less or equal than
$N(1-\pi)$. So, for the model to be consistent we require $\pi\leq v$.

We are interested in how the reproduction number $R_{0,e}=R_{0,e}(v,\pi)$
(where the $e$ stands for effective) changes once the vaccine is introduced in
the population and depends on $v$ and $\pi$. To analyse this, we compute the
dominant eigenvalue of the next generation matrix

$A^{\prime}=A^{\prime}(v,\pi)=\begin{bmatrix}\lambda_{G}(1-v)&(n_{c}-1)\frac{(1-v)}{(1-\pi)}\\\
\lambda_{G}(1-v)&0\end{bmatrix}.$

The results are depicted in Figure 1.

Figure 1: $R_{0,e}$ as a function of $v$, if both vaccines and mandates are
introduced.

In a class without a vaccine mandate, the distribution of the number $X$ of
unvaccinated children is $Bin(n_{c},\frac{1-v}{1-\pi})$ and we write
$\pi_{k}=P(X=k)$. In what follows we use

$u:=\frac{1-v}{1-\pi}.$

The expected number of unvaccinated children per class is given by
$\mu_{c}=n_{c}u=n_{c}\frac{1-v}{1-\pi}$. We observe that the probability that
a child, chosen uniformly at random among the unvaccinated children in the
population, belongs to a class with $k$ unvaccinated children is given by the
size-biased distribution, determined through

$\bar{\pi}_{k}=\frac{k\pi_{k}}{E[X]}=\frac{k\binom{n_{c}}{k}(u)^{k}(1-u)^{n_{c}-k}}{n_{c}u}=\binom{n_{c}-1}{k-1}u^{k-1}(1-u)^{n_{c}-1-k+1}.$
(3)

We notice that these are the probability masses of $Y+1$, where $Y=Y(v)$ is a
random variable with a binomial distribution with parameters $n_{c}-1$ and
$u$.

Using $\lambda_{G}$ from (2), we compute $R_{*,e}=R_{*,e}(v,\pi)$ which is the
household reproduction number after vaccination for the household model [3,
4], and represents the mean number of clumps (or in our case, classes) that
are infected by a typical infectious clump/class. It is defined to be
$R_{*,e}=\lambda_{G,v}\mu$, where $\lambda_{G,v}$ is the expected number of
global infections caused by a typical infected individual, and $\mu=\mu(v)$ is
the mean number of infected children in an infected class. $\lambda_{G,v}$ is
given by $\lambda_{G}(1-v)$, while

$\mu=\mathbb{E}[\bar{X}]=\sum_{k=1}^{n_{c}}k\bar{\pi}_{k}=1+(n_{c}-1)u.$

The latter follows from the size biased distribution in (3). Substituting this
in the expression for $R_{*,e}$, we obtain

$R_{*,e}=[1+(n_{c}-1)u](1-v)\lambda_{G}.$ (4)

The household reproduction number $R_{*,e}$, as well as most other
reproduction numbers defined in literature, including $R_{0,e}$, is useful
because an epidemic goes extinct with probability 1 if and only if
$R_{*,e}\leq 1$ [4]. The minimal value of $v$ for which $R_{*,e}\leq 1$, is
the critical vaccination coverage $v_{cr}$.

$R_{*,e}>1\Longleftrightarrow\frac{(1-v)^{2}}{(1-\pi)}\lambda_{G}(n_{c}-1)>1-(1-v)\lambda_{G}.$

Since $n_{c}$ represents the number of children in a class, the left hand side
is always non-negative. So, if the expected number of contacts with
unvaccinated children (i.e., $(1-v)\lambda_{G}$) is greater than $1$, i.e.
$1-v>\frac{1}{\lambda_{G}}$, then the epidemic is above threshold. Assume
instead that both terms are positive; then this inequality gives us an upper
bound for the fraction $1-\pi$ of schools without mandate:

$1-\pi<\frac{(1-v)^{2}\lambda_{G}(n_{c}-1)}{1-(1-v)\lambda_{G}}.$

So, if the number of schools without mandate is small enough and the
unvaccinated children are concentrated in those schools, the value of
$R_{*,e}$ is above 1 (see Figure 2).

Figure 2: Left: plot of $R_{*,e}$ as a function of $v\in[0.8,1]$ when
$R_{0}=15$, $n_{c}=25$ and $\pi=0$. The red asterisk denotes when $R_{*,e}=1$
and we call $v_{cr}$ the correspondent vaccination coverage.
Right: plot of $R_{*,e}$ as a function of $\pi$ when $v=v_{cr}$, again for
$R_{0}=15$ and $n_{c}=25$.

###  2.1 Branching process approximation

####  2.1.1 The branching process

Assume that the epidemic starts with one randomly chosen unvaccinated child in
the population and say that the class of the child has $\ell$ non-vaccinated
children. Then the probability that a class with $k$ unvaccinated children
gets infected by the class of the initially infected child is given by the
following expression

$\frac{\pi_{k}[1-(1-p_{G})^{k\ell}]}{\sum_{k^{\prime}=0}^{n_{c}}\pi_{k}^{\prime}[1-(1-p_{G})^{k^{\prime}\ell}]}.$
(5)

When $p_{G}<<1$ (as we assume), the expression $1-(1-p_{G})^{k\ell}$ is well
approximated by $k\ell p_{G}$, and we obtain that (5) is approximated by

$\frac{\pi_{k}p_{G}k\ell}{\sum_{k^{\prime}=0}^{n_{c}}\pi_{k}^{\prime}p_{G}k^{\prime}\ell}=\frac{\pi_{k}k}{\sum_{k^{\prime}=0}^{n_{c}}\pi_{k}^{\prime}k^{\prime}}=\bar{\pi}_{k}.$

If the number of classes is large, also the population size $N$ and the number
of unvaccinated children $N_{u}$ are large. So, we can use a branching process
approximation to describe the initial stages of the epidemic: in those early
stages infectious individuals contact only individuals that – with large
probability – are still susceptible, so the number of new infections caused by
households are asymptotically independent and global contacts lead to a tree-
structured epidemic generated graph, which has the law of a Galton-Watson tree
with offspring distribution described by $\bar{\pi}_{k}$
($k\in\mathbb{N}_{0}$) [6].

More formally, by a birthday problem approach [9, p.24], we can argue that the
probability that the first $k$ global contacts, $k=o(\sqrt{N_{u}})$, are all
with different households. Furthermore, it is well known (e.g. [1, Chapter 3])
that if we denote the number of infected children in an SIR epidemic by
$\lvert\mathcal{E}_{N_{u}}\rvert$ and the size of the approximating branching
process (i.e. a Galton-Watson process with offspring distribution described by
$\bar{\pi}_{k}$ ($k\in\mathbb{N}_{0}$)) by $\lvert\mathcal{B}\rvert$, then as
$N_{u}\to\infty$

$\mathbb{P}(\lvert\mathcal{E}_{N_{u}}\rvert>\log
N_{u})\to\mathbb{P}(\lvert\mathcal{B}\rvert=\infty).$

If $\lvert\mathcal{E}_{N_{u}}\rvert>\log_{u}N$ we talk about a major outbreak
of the epidemic.

####  2.1.2 The probability of avoiding infection for a medically exempt child

If the first infectious child belongs to a class of size $k$ (i.e. a class
with $k$ unvaccinated children), then the probability that the process dies
out is equal to the probability that the offspring of each student in the
first infectious class dies out. Since the infection rate for global contacts
between unvaccinated is $\lambda_{G}(1-v)$, a class of size $k$ will make
global contacts at a rate $\lambda_{G}(1-v)k$. This leads to the following
expressions for $\bar{q}$, the extinction probability of the approximating
branching process for an epidemic initiated by a uniformly chosen unvaccinated
individual from the population.

$\displaystyle\bar{q}$
$\displaystyle=\sum_{k=1}^{n_{c}}\mathbb{P}(\bar{X}=k)\sum_{\ell=0}^{\infty}\frac{[\lambda_{G}(1-v)k]^{\ell}}{\ell!}e^{-\lambda_{G}(1-v)k}\bar{q}^{l}$
(6)
$\displaystyle=\sum_{k=1}^{n_{c}}\mathbb{P}(\bar{X}=k)e^{-\lambda_{G}(1-v)k(1-\bar{q})}$
$\displaystyle=\sum_{j=0}^{n_{c}-1}\mathbb{P}(\bar{X}-1=j)e^{-\lambda_{G}(1-v)(j+1)(1-\bar{q})}$
$\displaystyle=\sum_{j=0}^{n_{c}-1}\binom{n_{c}-1}{j}u^{j}(1-u)^{n_{c}-1-j}e^{-\lambda_{G}(1-v)(j+1)(1-\bar{q})}$
$\displaystyle=(1-u+ue^{-\lambda_{G}(1-v)(1-\bar{q})})^{n_{c}-1}e^{-\lambda_{G}(1-v)(1-\bar{q})}$

which means that $\bar{q}$ is a solution of the fixed point equation

$\bar{q}=(1-u+ue^{-\lambda_{G}(1-v)(1-\bar{q})})^{n_{c}-1}e^{-\lambda_{G}(1-v)(1-\bar{q})}.$

From the theory of branching processes [11], we know that $\bar{q}$ is the
smallest solution of this equation, which is strictly less than 1 if and only
if $R_{*,e}>1$.

By assumption a medically exempt child is enrolled in a school that has a
vaccine mandate. This means that they belong to a class of effective size 1.
By the epidemic generated graph representation of an epidemic [6], we know
that $p_{me}$, the probability that a medically exempt child avoids being
infected if a large outbreak occurs, is equal to the probability that a large
outbreak occurs if the outbreak starts with a single infected child in a class
of effective size 1. This result can also be obtained by using susceptibility
sets [2] and see that the so-called forward and backward branching process
have the same offspring distribution. Using arguments as in (6), we obtain
that

$p_{me}=p_{me}(v,\pi)=e^{-\lambda_{G}(1-v)(1-\bar{q})}.$ (7)

If there is no vaccine mandate, a medically exempt child is just like any
other unvaccinated child and the probability of avoiding infection is
$\bar{q}$ for the case that $\pi=0$. We denote this probability by
$p_{wm}=p_{wm}(v)$, where wm stands for “without mandate”.

##  3 Results

Considering the entire population, we may expect that mandates make the global
situation worse, because unvaccinated children have on average more other
unvaccinated children in their class and children in the same class contact
each other at a higher rate than other pairs of children. But in this
situation a medically exempt child can only be infected by global contacts. In
the case without mandate, even if the global situation is more controlled,
such a child could be enrolled in a class with other unvaccinated children.
And the probability $p_{me}$ for a given medically exempt child to avoid
infection in a situation with vaccine mandates may be larger than that
probability in a situation without mandates $p_{wm}$. In this section we
analyse the model to justify this intuition.

In Figure 3 we consider the case where $\pi=1/2$ and $v$ is varying. Here
$p_{me}$ is higher than $p_{wm}$ up to $v\leq v_{s}\approx 93\%$, and we have
the opposite inequality for $v>v_{s}$, where the $s$ stands for switching
level. The graph tells us that if $v$ is low, a medically exempt child is
profiting from vaccine mandates at schools, since they will for sure be
enrolled in a class with no other unvaccinated children. But when the vaccine
percentage is sufficiently high the probability of having only one
unvaccinated child per class is already high and mandates worsen the global
situation for everyone. Note that $p_{wm}$ reaches $1$ at a lower $v$ value
than $p_{me}$. This confirms the intuition that without mandates a lower
vaccine coverage is sufficient to guarantee herd-immunity.

Figure 3: Plot of the probability $p_{me}$ that a medically exempt child
avoids infection as a function of $v$ with $\pi$ set to $0.5$, and of the
probability $p_{wm}$ that a random unvaccinated child avoids infection in the
case with no mandates ($\pi=0$). We see that the two lines cross at
$v_{s}\approx 0.925$. The minimal value of $v$ for which $p_{wm}=1$ is at
$v\approx 0.933$, while the minimal value of $v$ for which $p_{me}=1$ is at
$v\approx 0.95$.

In Figure 4, we again consider $p_{me}$, the probability that a medically
exempt individual avoids infection, but now as a function of $\pi$, when
$v=0.93$ and $v=0.90$. For $v=0.93$, we can see that, for low values of $\pi$
(lower than approximately $0.20$), the medically exempt are more protected in
the model with mandates. But if $\pi$ increases, then the mandates are
actually disadvantageous to the medically exempt children, even if they are
the only unvaccinated children in their class. If $v=0.90$, already without
any vaccine mandates, the global spread of the disease is so easy, that the
global spread does not increase much by introducing vaccine mandates, while
the relative protection of having no other unvaccinated children in the class
is considerable and vaccine mandates are always beneficial for medically
exempt children.

Figure 4: Plot of $p_{me}$, the probability that a medically exempt child
avoids infection as a function of $\pi$, when $v=0.93$ (left) and when
$v=0.90$ (right). We compare it with the probability $p_{wm}$ that a random
unvaccinated child avoids infection in the case with no mandates ($\pi=0$).

Up to now, we have assumed that the introduction of mandates in a society will
not change the vaccine percentage. This is not realistic, since some parents
might vaccinate their children in order to let them go to the school of their
preference. We call those parents “converted”. We may think of $v$ as the sum
of two components: $v_{i}$ that is the initial proportion of vaccinated
individuals and $v_{c}$ that is the additional portion of “converted”. Note
that $v_{c}\leq 1-v_{i}$. The graph in Figure 5 shows a plot of $R_{*,e}$ as a
function of $v_{c}$, when $\pi=0.5$, $v_{i}=0.9$ and $n_{c}=25$.

Figure 5: Left: $R_{*,e}$ as a function of $v_{c}$. The red asterisk is the
intersection point with $y=1$.
Right: Plot of $p_{me}$ and $p_{wm}$ as a function of $v_{c}$. Note that this
figure may be deduced from Figure 3 with $v$ restricted to $[0.9,1]$. This is
because we fix the initial fraction of vaccinated $v_{i}$ and then define
$v=v_{i}+v_{c}$ and compute $p_{me}$ and $p_{wm}$ with such combined $v$.

##  4 Discussion

In this paper we study how vaccine mandates might affect the spread of an
infectious disease in elementary schools, and especially how the probability
of being infected changes for medically exempt children. This effect may be
counter-intuitive and local mandates might worsen the risks for medically
exempt children. This work was inspired by a law that was proposed by the
Dutch House of Representatives and was later rejected by the Senate [8]. In
this measure, childcare owners could freely decide whether or not joining the
public vaccination program was compulsory for children to enrol at the
daycare.

In our study we do not distinguish between child-daycare and elementary
schools, and we only allow the disease to spread within schools, while the
population structure and contact patterns are actually much more complex and
the infection could take place in households, sport centers and so on.
Moreover, we do not consider multiple waves of infection nor the possibility
of becoming susceptible after being recovered, which indeed is not that
relevant for measles, but might be for other diseases.

Despite the limitations due to the model hypotheses, the analysis describes
the qualitative different behaviours of the two situations (with and without
vaccine mandates) and can be important for further ethical and political
decisions made on vaccines.

The analysis can easily be extended to diseases with lower reproduction number
than measles, see Appendix C. It depends on the relative spread within and
between households, whether this will lead to less or more pronounced possible
undesirable effects of vaccine mandates.

##  Appendix A Random class sizes

Our analysis is easily extended to random class sizes, in which the number of
children in a class is distributed as a random variable $N_{c}$. Call
$p_{k}=\mathbb{P}(N_{c}=k)$; then the probability that a newly infected class
is of size $k$ is given by the size-biased version $\bar{N_{c}}$ of $N_{c}$
and we define $\bar{p}_{k}=\mathbb{P}(\bar{N_{c}}=k)$. To evaluate the
expected number of type-$2$ individuals that a type-$1$ individual infects
(which is element $(1,2)$ of matrix $A$ in (1)), we compute

$\mathbb{E}[\bar{N_{c}}-1]=\sum_{k=1}^{\infty}(k-1)\bar{p}_{k}=\sum_{k=1}^{\infty}\frac{k(k-1)p_{k}}{\mathbb{E}[N_{c}]}=\frac{\mathbb{E}[N_{c}(N_{c}-1)]}{\mathbb{E}[N_{c}]}.$

For example, if $N_{c}\sim Poiss(m)$, say with $m=25$. Then
$\mathbb{E}[N_{c}]=m$ and, with the same procedure as in Section 2, we obtain
$\lambda_{G}\approx 5.6$.

We denote the random number of unvaccinated per class after both vaccines and
mandates are introduced in the model again by $X$. Then $X$ conditioned on
$N_{c}=k$ is $Bin(k,u)$ distributed and

$\mathbb{P}(X=l)=\sum_{k}\mathbb{P}(X=l|N_{c}=k)\mathbb{P}(N_{c}=k)$

while

$\mathbb{P}(\bar{X}=l)=\frac{l\mathbb{P}(X=l)}{\mathbb{E}[X]}=\sum_{k}\mathbb{P}(\bar{X}=l|N_{c}=k)\mathbb{P}(N_{c}=k).$
(8)

We observe that

$\mathbb{P}(\bar{X}=l|N_{c}=k)=\frac{l\mathbb{P}(X=l|N_{c}=k)}{\mathbb{E}[X|N_{c}=k]}$

which implies, after some computations, that $\bar{X}$ conditioned on
$N_{c}=k$ is distributed as $1+Bin(k-1,u)$, and then allows us to compute the
unconditioned distribution of $\bar{X}$ using (8).

##  Appendix B Non trivial local infection probabilities

Throughout all our analysis above, we have considered the case in which the
local infection probability $p_{L}=1$, which simplifies the model as explained
in Section 1.2. Here we want to study non trivial cases with $p_{L}<1$, that
correspond to a different contact pattern.

Assume that an infectious individual infects their classmates with probability
$p_{L}=1/2$ per classmate during their infectious period. Setting both $v$ and
$\pi$ to zero, we compute the global infection rate $\lambda_{G}$ as in
Section 2. We consider a class of size $25$ in which one student gets infected
from a global contact. The probability that a given one of their classmates
avoids infection within two steps is
$\frac{1}{2}\Big{(}1-\frac{1}{4}\Big{)}^{23}$ and the probability that at
least one of the classmates avoids infection is at least
$24\times\frac{1}{2}\Big{(}1-\frac{1}{4}\Big{)}^{23}\approx 0$. This means
that all the pupils of that class will have contracted the virus within two
generations of infection with very high probability. A primary case infects on
average $(n_{c}-1)/2$ classmates which will, in turn, infect as many remaining
susceptibles. The next generation matrix (1) then becomes

$A=\begin{bmatrix}\lambda_{G}&\frac{n_{c}-1}{2}&0\\\ \lambda_{G}&0&1\\\
\lambda_{G}&0&0\\\ \end{bmatrix}$

and $R_{0}$ solves

$R_{0}^{2}(\lambda_{G}-R_{0})+R_{0}\lambda_{G}\frac{n_{c}-1}{2}+\lambda_{G}\frac{n_{c}-1}{2}=0.$

If $R_{0}=15$ and $n_{c}=25$ then

$\lambda_{G}=\frac{R_{0}^{3}}{R_{0}^{2}+\frac{n_{c}-1}{2}+R_{0}\frac{n_{c}-1}{2}}\approx
8.1.$

Recall that if $p_{L}=1$, then $\lambda_{G}\approx 5.8$.

Now introduce both vaccines and mandates in the model. We consider a class
with $n+1$ unvaccinated children, one of which is infected through a global
contact. The epidemic that starts spreading via local contacts in the class is
described by – to use the notation of Andersson and Britton [1, Chapter 2] –
the standard SIR epidemic model $E_{n,1}(\lambda_{L},1)$, where $\lambda_{L}$
indicates the local infectious rate. To compute the value of $\lambda_{L}$ we
link it to $p_{L}$ by observing that an infectious individual makes contacts
with a given classmate at the time steps of a
Poisson$\big{(}\lambda_{L}/n\big{)}$ process; therefore, the probability
$1-p_{L}$ of having no contacts with a given individual is given by
$e^{-\lambda_{L}/n}$.

We indicate by $P^{n}_{j}$ the probability that $j$ of the $n$ initial
susceptibles get infected during the (local) epidemic. The distribution of
these probabilities is obtained by solving a lower triangular linear system
given in [1, Chapter 2, Theorem 2.2]. Then, in this case, the probability of
extinction $\bar{q}$ in (6) will be

$\displaystyle\bar{q}$
$\displaystyle=\sum_{k=1}^{n_{c}}\mathbb{P}(\bar{X}=k)\sum_{j=0}^{k-1}P^{k-1}_{j}\sum_{l=0}^{\infty}\frac{[\lambda_{G}(1-v)(j+1)]^{l}}{l!}e^{-\lambda_{G}(1-v)(j+1)}\bar{q}^{l}$
(9)
$\displaystyle=\sum_{k=1}^{n_{c}}\mathbb{P}(\bar{X}=k)\sum_{j=0}^{k-1}P^{k-1}_{j}e^{-\lambda_{G}(1-v)(j+1)(1-\bar{q})}.$

We observe that, keeping $p_{L}$ fixed and increasing the number of
susceptibles (of the local epidemic) $n$, the probability $P^{n}_{n}$ that
everyone gets infected gets closer and closer to $1$, as seen in Figure 6.

In Figure 7 we see, on the left, a plot of the probabilities $p_{me}$ and
$p_{wm}$ that a medically exempt child avoids infection in the case with and
without class mandates respectively, when $\pi=1/2$, as functions of $v$. On
the right, we can see a comparison between the cases $p_{L}=1$ and
$p_{L}=1/2$. We observe that, even if in the latter case the local contact
probability is halved, the probability $p_{me}$ that a medically exempt child
avoids infection is higher in the first case. This depends on the value of the
global infection rate $\lambda_{G}$, which we remember increases from 5.8 in
the first case to 8.1 in the second: if $R_{0}$ stays the same but $p_{L}$
decreases, then $p_{G}$ increases, which means that the medically exempts are
more exposed to the virus. For the probability $p_{wm}$ that a random
unvaccinated avoids infection in the case without mandates we conversely see
almost no difference in the two cases.

Figure 6: Value of $P^{n}_{n}$ plotted as a function of $n$.

Figure 7: Left: Plot of the probability $p_{me}$ that a medically exempt child
avoids infection as a function of $v$ with $\pi$ set to $0.5$, and of the
probability $p_{wm}$ that a random unvaccinated child avoids infection in the
case with no mandates ($\pi=0$), when the local probability $p_{L}$ is set to
$1/2$. The two lines cross at $v_{s}\approx 0.925$, as in Figure 3.

Right: Confront between the cases $p_{L}=1$ (see Figure 3) and $p_{L}=1/2$.

Similarly we can analyse the cases $p_{L}=1/3$ and $p_{L}=1/5$. For the same
reason as above, we can assume each individual in a class to be infected
within a maximum of two generations of infectives. If one child is infected by
a global contact, it will infect on average $(n_{c}-1)/3$ and $(n_{c}-1)/5$
classmates respectively, which will in turn infect the remaining
$2(n_{c}-1)/3$ and $4(n_{c}-1)/5$ children. The next generation matrices
become

$A_{1/3}=\begin{bmatrix}\lambda_{G}&\frac{n_{c}-1}{3}&0\\\ \lambda_{G}&0&2\\\
\lambda_{G}&0&0\\\ \end{bmatrix}\qquad
A_{1/5}=\begin{bmatrix}\lambda_{G}&\frac{n_{c}-1}{5}&0\\\ \lambda_{G}&0&4\\\
\lambda_{G}&0&0\\\ \end{bmatrix}$

which lead to $\lambda_{G}\approx 9.3$ and 10.7 respectively. Then, using the
same expression for $\bar{q}$ found in (9), we obtain similar graphs for
$p_{me}$ and $p_{wm}$, as seen in Figure 8.

Finally, in Figure 9 we can see a comparison between the plots for $p_{me}$
(on the left) and $p_{wm}$ (on the right) for the different $p_{L}$ values we
have previously considered. What we observe is that, while $p_{me}$ tends to
decrease when $p_{L}$ decreases, $p_{wm}$ has the opposite behaviour. When
mandates are not present and unvaccinated are uniformly spread in the
classrooms, having a lower $p_{L}$ value slightly increases the probability
$p_{wm}$ of avoiding infection, as the clusters of infectives forming during
the spread of the virus are smaller. When mandates are present, the
probability $p_{me}$ that a medically exempt child avoids infection decreases
with $p_{L}$, because the global infection rate $\lambda_{G}$ reaches higher
values.

Figure 8: Plot of the probability $p_{me}$ that a medically exempt child
avoids infection as a function of $v$ with $\pi$ set to $0.5$, and of the
probability $p_{wm}$ that a random unvaccinated child avoids infection in the
case with no mandates ($\pi=0$), when the local probability $p_{L}$ is set to
$1/2$ (left) and $1/5$ (right).

Figure 9: Comparison of the plots of the probability $p_{me}$ that a medically
exempt child avoids infection as a function of $v$ when $\pi=0.5$ (left), and
of the probability $p_{wm}$ that a random unvaccinated child avoids infection
in the case of no vaccine mandates (right), for the different $p_{L}$ values
considered above.

##  Appendix C Less infectious diseases

We are now interested in considering the case of a less infectious disease;
for example, we take the basic reproduction number $R_{0}=6$.

Comparing Figure 3 and Figure 10 we see that in the first case the functions
describing $p_{me}$ and $p_{wm}$ increase much faster than in the case
$R_{0}=6$. This means that, with a less infectious virus, the difference
between the cases with and without mandates is more pronounced. In fact, as
shown in Figure 10, for a medically exempt child escaping infection is
significantly more probable in the case with vaccine mandates, up to the
switching point $v\approx 0.82$.

Figure 10: Plot of the probability $p_{me}$ that a medically exempt child
avoids infection as a function of $v$ with $\pi$ set to $0.5$, and of the
probability $p_{wm}$ that a random unvaccinated child avoids infection in the
case with no mandates ($\pi=0$). We see that the two lines cross at
$v_{s}\approx 0.82$. The minimal value of $v$ for which $p_{wm}=1$ is at
$v\approx 0.83$, while the minimal value of $v$ for which $p_{me}=1$ is at
$v\approx 0.88$. To compare it with the case $R_{0}=15$ see Figure 3.

In Figure 11 we see how $p_{me}$ and $p_{wm}$ vary as functions of $\pi$ when
the value of $v$ is fixed and close to the herd immunity; the behaviour is the
same as the one observed in Figure 4.

Figure 11: Plot of $p_{me}$, the probability that a medically exempt child
avoids infection as a function of $\pi$, when $v=0.83$ (left) and when
$v=0.80$ (right). We compare it with the probability $p_{wm}$ that a random
unvaccinated child avoids infection in the case with no mandates ($\pi=0$).
See Figure 4 to compare it with the case $R_{0}=15$.

## References

  * [1] H. Andersson and T. Britton.  Stochastic epidemic models and their statistical analysis, volume 151.  Springer Science & Business Media, 2000.
  * [2] F. Ball and O. D. Lyne.  Stochastic multi-type sir epidemics among a population partitioned into households.  Advances in Applied Probability, 33(1):99–123, 2001.
  * [3] F. Ball, D. Mollison, and G. Scalia-Tomba.  Epidemics with two levels of mixing.  The Annals of Applied Probability, 7(1):46–89, 1997.
  * [4] F. Ball, L. Pellis, and P. Trapman.  Reproduction numbers for epidemic models with households and other social structures II: Comparisons and implications for vaccination.  Mathematical Biosciences, 274:108–139, 2016.
  * [5] N. G. Becker and K. Dietz.  The effect of household distribution on transmission and control of highly infectious diseases.  Mathematical Biosciences, 127(2):207–219, 1995.
  * [6] T. Britton, K. Y. Leung, and P. Trapman.  Who is the infector? general multi-type epidemics and real-time susceptibility processes.  Advances in Applied Probability, 51(2):606–631, 2019.
  * [7] O. Diekmann, H. Heesterbeek, and T. Britton.  Mathematical tools for understanding infectious disease dynamics, volume 7.  Princeton University Press, 2013.
  * [8] Eerste Kamer (Dutch Senate).  Initiatiefvoorstel-Raemakers en Van Meenen, 2022.
  * [9] G. Grimmett and D. Stirzaker.  Probability and random processes.  Oxford university press, third edition, 2002.
  * [10] J. Hübschen, I. Gouandjika-Vasilache, and J. Dina.  Measles.  Lancet, 399:678–690, 2022.
  * [11] P. Jagers.  Branching processes with biological applications.  Wiley, 1975.
  * [12] K. Kuulasmaa.  The spatial general epidemic and locally dependent random graphs.  Journal of Applied Probability, 19(4):745–758, 1982.
  * [13] Ministero della Salute (Italian Ministery of Health).  Vaccinazioni, 2021.
  * [14] L. Pellis, F. Ball, and P. Trapman.  Reproduction numbers for epidemic models with households and other social structures. I. definition and calculation of ${R}_{0}$.  Mathematical Biosciences, 235(1):85–97, 2012.
  * [15] S. Thompson, J. C. Meyer, R. J. Burnett, and S. M. Campbell.  Mitigating vaccine hesitancy and building trust to prevent future measles outbreaks in england.  Vaccines, 11(2):288, 2023.
  * [16] World health Organisation.  History of measles vaccine.  https://www.who.int/news-room/spotlight/history-of-vaccination/history-of-measles-vaccination.
  * [17] World Health Organisation.  Measles.  https://www.who.int/news-room/fact-sheets/detail/measles, March 2023\.


# A note on extremal Sombor indices
of trees with a given degree sequence

Ivan Damnjanović
University of Niš, Faculty of Electronic Engineering, Niš, Serbia
<EMAIL_ADDRESS>
Diffine LLC, San Diego, California, USA
<EMAIL_ADDRESS>ID is supported by Diffine LLC. Marko Milošević
University of Niš, Faculty of Sciences and Mathematics, Niš, Serbia
<EMAIL_ADDRESS>Dragan Stevanović
Mathematical Institute of the Serbian Academy of Sciences and Arts,
Belgrade, Serbia
<EMAIL_ADDRESS>DS is supported by the Serbian Ministry of
Education, Science and Technological Development through the Mathematical
Institute of SASA, and by the project F-159 of the Serbian Academy of Sciences
and Arts.

###### Abstract

We note here that the problem of determining extremal values of Sombor index
for trees with a given degree sequence fits within the framework of results by
Hua Wang from [Cent. Eur. J. Math. 12 (2014) 1656–1663], implying that the
greedy tree has the minimum Sombor index, while an alternating greedy tree has
the maximum Sombor index.

In a recent private communication, Ivan Gutman asked a number of colleagues to
characterize tree(s) of order $n$ with the given degree sequence $\mathcal{D}$
whose Sombor indices are minimum and maximum. Recall that Sombor index $SO(G)$
of a graph $G=(V,E)$ is defined in [2] as

$SO(G)=\sum_{uv\in E}\sqrt{d_{u}^{2}+d_{v}^{2}},$

where $d_{u},d_{v}$ are degrees of the vertices $u,v\in V$. Observing that
switching edges (i.e., delete edges $ab$ and $cd$, and add edges $ac$ and
$bd$) decreases Sombor index under suitable conditions, while keeping degrees
intact, we quickly jumped in to show in [1] that a greedy tree necessarily
attains the minimum Sombor index among trees with degree sequence
$\mathcal{D}$. Actually, the greedy tree is the unique tree that minimizes
pseudo-Sombor index (see [1] for details), but in principle there may exist
other trees with the same minimum value of Sombor index as the greedy tree for
given $\mathcal{D}$. A more detailed reading of references by one of us during
the subsequent attempt to determine trees with the maximum Sombor index,
revealed that this problem actually fits within the framework of results by
Hua Wang [3], which quickly implies both that the greedy tree has the minimum
and that an alternating greedy tree has the maximum value of Sombor index
among trees with degree sequence $\mathcal{D}$.

For the sake of completeness, let us present here the result of Wang [3].
Assume that the degree sequence $\mathcal{D}$ is ordered in a non-increasing
order, and denote by $d_{1}\geq\dots\geq d_{m}$ the degrees of internal
vertices (i.e., the elements of $\mathcal{D}$ that are greater than one). Both
the greedy tree and the alternating greedy tree are constructed
algorithmically. The greedy tree is constructed as follows:

  1. (g1)

Label the root with the largest degree $d_{1}$;

  2. (g2)

Label the neighbors of the root as $d_{2}$, $d_{3}$, …, assigning to each next
neighbor the largest available degree;

  3. (g3)

For each labelled vertex in the current level, considered in a non-increasing
order of labels, label its children in turn with the largest available degree;

  4. (g4)

Repeat (g3) as long as there are available internal degrees, then add
necessary number of leaves so that the degree of each vertex is equal to its
label.

The alternating greedy tree is constructed by a recursive procedure:

  1. (a1)

If $m-1\leq d_{m}$, the alternating greedy tree is a tree rooted at the vertex
$r$ with $d_{m}$ children, among which $d_{m}-m+1$ are leaves, while the
remaining children have degrees $d_{1},\dots,d_{m-1}$;

  2. (a2)

If $m-1\geq d_{m}+1$, create a subtree $T$ rooted at $r$ with $d_{m}-1$
children with degrees $d_{1},\dots,d_{d_{m}-1}$;

  3. (a3)

Let $S$ be the alternating greedy tree corresponding to the sequence
$(d_{d_{m}},\dots,\linebreak d_{m-1})$ of internal degrees, and let $v$ be a
leaf of $S$ with the smallest degree of its neighbor. The alternating greedy
tree for the sequence $(d_{1},\dots,d_{m})$ is obtained by identifying the
root $r$ of $T$ with $v$ in $S$.

Figure 1: Illustration of constructions of: (a) the greedy tree, (b) the
constituents of alternating greedy trees, and (c–e) some feasible alternating
greedy trees for the sequence $(5,4,3,3,3,2,2,2)$ of internal vertex degrees.

These two constructions are illustrated in Fig. 1 for the internal degree
sequence $(5,4,3,3,3,2,2,2)$. The greedy tree, shown in Fig. 1(a), is produced
uniquely by the steps (g1)–(g4). However, several non-isomorphic alternating
greedy trees may be produced by the steps (a1)–(a3), as it can happen that the
leaf $v$ in the step (a3) can be selected in different non-isomorphic ways.
Namely, step (a2) applied to $(5,4,3,3,3,2,2,2)$ produces the subtree $T_{1}$
with the root $r_{1}$, leaving the subsequence $(4,3,3,3,2,2)$ for which one
has to produce an alternating greedy tree in step (a3). This calls step (a2)
recursively to produce the subtree $T_{2}$ with the root $r_{2}$, and leaves
the subsequence $(3,3,3,2)$. Another recursive call to step (a2) produces the
subtree $T_{3}$ with the root $r_{3}$, which leaves the subsequence $(3,3)$
for which step (a1) produces the alternating greedy tree $S_{3}$. All these
“constituents” are shown in Fig. 1(b). However, going back from these
recursive calls and continuing with step (a3) yields several possible choices
for the choice of the leaf $v$. First, the root $r_{3}$ of $T_{3}$ may be
identified with either of the leaves $v_{1}$ and $v_{2}$ of $S_{3}$. If
$r_{3}$ is identified with $v_{1}$, as done in Fig. 1(c) and 1(d), then the
root $r_{2}$ of $T_{2}$ may be further identified with either one of the
remaining leaves of $S_{3}$ or one of the leaves of $T_{3}$ in the newly
formed tree. After this is done, there are still several choices left for the
choice of the leaf which should be identified with the root $r_{1}$ of
$T_{1}$. Fig. 1(c)–(e) shows some of the final resulting alternating greedy
trees (where in Fig. 1(e) the root $r_{3}$ was initially identified with the
leaf $v_{2}$).

Generalizing earlier results, Wang [3] considered for a tree $T=(V,E)$ the
general form of a topological index defined as

$R_{f}(T)=\sum_{uv\in E}f(d_{u},d_{v}),$

where $f\colon\mathbb{N}\times\mathbb{N}\to\mathbb{R}$ is a symmetric
function. He proved the following theorem.

######  Theorem 1 ([3]).

If the symmetric function $f\colon\mathbb{N}\times\mathbb{N}\to\mathbb{R}$
satisfies

$f(x,a)+f(y,b)\geq f(y,a)+f(x,b)\quad\mbox{for all $x\geq y$ and $a\geq b$}$
(1)

(with strict inequality implied if both $x>y$ and $a>b$), then $R_{f}(T)$ is
maximized by the greedy tree and minimized by an alternating greedy tree among
trees with given degree sequence.

Let us now specifically define $f$ as

$f(x,a)=-\sqrt{x^{2}+a^{2}}$

so that $R_{f}(T)$ is actually minus Sombor index of $T$. In this case the
condition (1) reads as

$-\sqrt{x^{2}+a^{2}}-\sqrt{y^{2}+b^{2}}\geq-\sqrt{y^{2}+a^{2}}-\sqrt{x^{2}+b^{2}}$

whenever $x\geq y$ and $a\geq b$, which is equivalent to

$\sqrt{x^{2}+a^{2}}+\sqrt{y^{2}+b^{2}}\leq\sqrt{y^{2}+a^{2}}+\sqrt{x^{2}+b^{2}}.$

After squaring and rearranging, this is further equivalent to

$(x^{2}+a^{2})(y^{2}+b^{2})\leq(y^{2}+a^{2})(x^{2}+b^{2}),$

and further to

$0\leq(a^{2}-b^{2})(x^{2}-y^{2}),$

which is certainly satisfied whenever $x\geq y$ and $a\geq b$ (with strict
inequality if both $x>y$ and $a>b$). Hence Theorem 1 holds for minus Sombor
index, leading to the following corollary for Sombor index itself.

######  Corollary 2.

Sombor index is minimized by the greedy tree and maximized by an alternating
greedy tree among trees with given degree sequence.

## References

  * [1] I. Damnjanović, D. Stevanović, Greedy trees have minimum Sombor indices, arXiv:2211.05559 (2022).
  * [2] I. Gutman, Geometric approach to degree-based topological indices: Sombor indices, MATCH Commun. Math. Comput. Chem. 86 (2021) 11–16.
  * [3] H. Wang, Functions on adjacent vertex degrees of trees with given degree sequence, Cent. Eur. J. Math. 12 (2014) 1656–1663.


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

M.T. Gallagher1,2,3<EMAIL_ADDRESS>
2School of Mathematics, 3Institute for Metabolism and Systems Research, and
4Institute of Translational Medicine, University of Birmingham, B15 2TT, UK
5John Tyndall Institute, School of Engineering, and 6School of Medicine and
Dentistry, University of Central Lancashire, Preston, PR1 2HE, UK R.A.J.
Wain4,5,6<EMAIL_ADDRESS>
2School of Mathematics, 3Institute for Metabolism and Systems Research, and
4Institute of Translational Medicine, University of Birmingham, B15 2TT, UK
5John Tyndall Institute, School of Engineering, and 6School of Medicine and
Dentistry, University of Central Lancashire, Preston, PR1 2HE, UK S. Dari2
J.P. Whitty5<EMAIL_ADDRESS>
2School of Mathematics, 3Institute for Metabolism and Systems Research, and
4Institute of Translational Medicine, University of Birmingham, B15 2TT, UK
5John Tyndall Institute, School of Engineering, and 6School of Medicine and
Dentistry, University of Central Lancashire, Preston, PR1 2HE, UK D.J.
Smith2,3<EMAIL_ADDRESS>
2School of Mathematics, 3Institute for Metabolism and Systems Research, and
4Institute of Translational Medicine, University of Birmingham, B15 2TT, UK
5John Tyndall Institute, School of Engineering, and 6School of Medicine and
Dentistry, University of Central Lancashire, Preston, PR1 2HE, UK

###### Abstract

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.

##  1 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:

$\mu(\dot{\gamma})=\mu_{\infty}+\frac{\mu_{0}-\mu_{\infty}}{1+(\lambda\dot{\gamma})^{1-n}},$
(1)

where $\mu$ is the effective viscosity of the fluid as a function of shear
rate $\dot{\gamma}$, the parameters $\mu_{0}$ and $\mu_{\infty}$ are the zero
and infinite limit shear viscosities respectively, $\lambda$ 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,

$\mu\left(\dot{\gamma}\right)=\frac{\mu_{0}}{\left(1+\left(\lambda\dot{\gamma}\right)^{2}\right)^{\frac{1-n}{2}}}.$
(2)

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:

$\mu\left(\dot{\gamma}\right)=\frac{\mu_{0}}{\left(1+\left(\lambda\dot{\gamma}\right)^{a}\right)^{\frac{1-n}{a}}}.$
(3)

This model has four free parameters ($\mu_{0},\lambda,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

$\mu\left(\dot{\gamma}\right)=\mu_{\infty}+\frac{\mu_{\infty}-\mu_{0}}{\left(1+\left(\lambda\dot{\gamma}\right)^{a}\right)^{\frac{1-n}{a}}}.$
(4)

Equation (4) differs from both Carreau’s and Yasuda’s models through the
inclusion of an infinite shear rate viscosity parameter $\mu_{\infty}$ (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.

The variability in terminology can also be found in major commercial
computational fluid dynamics codes such as ANSYS-Fluent, ANSYS-CFX, and Abaqus
(see Table 1). For the avoidance of confusion, below we refer to equations
(1), (2), (3) and (4) as the Cross-1965, Carreau-1972, Yasuda-1979 and
BCCY-1987 (Bird-Cross-Carreau-Yasuda) models respectively.

Table 1: Three commercial CFD packages and their varying terminology Software | Their name | Equation
---|---|---
Abaqus | Carreau-Yasuda |  Eq (4)
| Carreau |  Eq (4) with $a=2$
| Cross |  Eq (1)
ANSYS-CFX | Carreau-Yasuda |  Eq (4)
| Bird-Carreau |  Eq (4) with $a=2$
| Cross |  Eq (8)
ANSYS-Fluent | Carreau |  Eq (4) with $a=2$ and a
|  | leading temperature factor
| Cross |  Eq (8) with a leading
|  | temperature factor

##  2 Data fitting

Given a functional form $\mu(\dot{\gamma};\bm{\theta})$ for the viscosity
$\mu$ at shear rate $(\dot{\gamma})$ with parameters
$\bm{\theta}=(\mu_{0},n,\ldots)$, and rheometry data
$(\dot{\gamma}_{m},\mu_{m})$, maximum likelihood estimation with the
assumption of normally distributed error leads to the least squares parameter
estimate,

$\bm{\theta}^{*}=\underset{\bm{\theta}}{\mathrm{argmin}}\
\mathcal{L}(\bm{\theta}):=\underset{\bm{\theta}}{\mathrm{argmin}}\sum_{m=1}^{M}(\mu(\dot{\gamma}_{m};\bm{\theta})-\mu_{m})^{2}.$
(5)

We will also find it useful to optimize over subsets of parameters while the
remaining are held fixed. To denote this we will use, for example, the
notation

$\mathcal{L}(\mu_{0},\lambda,*):=\underset{n}{\mathrm{max}}\
\mathcal{L}(\mu_{0},\lambda,n).$ (6)

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,

$\bm{\theta}^{\Diamond}=\underset{\bm{\theta}}{\mathrm{argmin}}\
\mathcal{\ell}(\bm{\theta}):=\underset{\bm{\theta}}{\mathrm{argmin}}\sum_{m=1}^{M}(\log(\mu(\dot{\gamma}_{m};\bm{\theta}))-\log(\mu_{m}))^{2}.$
(7)

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.

##  3 Results

Each model will be considered in turn, starting with the model possessing the
fewest free parameters. The units of $\mu_{0},\ \mu_{\infty}$ will be
centipoise, $\lambda$ will be seconds, and $n,a$ are dimensionless.

###  3.1 Carreau-1972 three-parameter model

Performing a maximum likelihood fit of the Carreau-1972 model with normal
errors, i.e. by minimizing $\mathcal{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 $\mu_{0}$ or the maximum value of $\lambda$ depending on
the limits chosen, however with $n=0.483$ consistently. The reason for this
behaviour is evident in Figure 1; the likelihood surface
$\mathcal{L}(\mu_{0},\lambda,*)$ exhibits an extended ‘ridge’. Taking an upper
bound of $\mu_{0}<350$ yields the parameter estimate $\mu_{0}=301,\
\lambda=200,\ n=0.483$, depicted with a blue star; the cost function value is
$\mathcal{L}^{*}=25.8$. To show how indeterminate this fit is, we will examine
an arbitrarily chosen ‘alternative’ parameter tuple of $\mu_{0}=211,\
\lambda=100$ and fitted optimum $n^{+}=0.482$ (red plus), which has a very
similar cost function value of $\mathcal{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 $\mu_{0}$ or $\lambda$. It is also
of note that either parameter set fits the data well with these parameters up
to approximately $\dot{\gamma}=10$ s-1, however they both perform rather badly
for higher values of shear rate.

Figure 1: Carreau-1972 model fits (normal error) to data of Skalak et al.
Panel (a) shows linearly spaced contours of the objective function
$\mathcal{L}\left(\mu_{0},\lambda,*\right)$, with a blue star showing the
location of the Matlab fmincon parameter fit (with $n=0.483$), and a red plus
showing the arbitrarily chosen ‘alternative’ parameter choice (with
$n^{+}=0.482$). Panels (b) and (c) plot the fits of the Carreau-1972 model to
the data of Skalak et al. (1981) for these parameter tuples shown on linear
and logarithmic axes respectively.

One may ask whether the more traditional approach of fitting to the log-log
plot, i.e. by taking lognormal error and minimizing $\ell$, 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 $\mu_{0}^{\Diamond}=137,\ \lambda^{\Diamond}=200,\ n^{\Diamond}=0.635$
which has cost function $\ell^{\Diamond}=0.732$. A manually and arbitrarily-
chosen tuple $\mu_{0}^{\times}=107,\ \lambda^{\times}=100,\ n^{\times}=0.634$,
plotted as a red cross, yields a very similar cost function value of
$\ell^{\times}=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 $\mu_{0},\lambda$ are again indeterminate
from the data. The same issue occurs for several other experimental blood
rheology data sets, see A).

Figure 2: Carreau-1972 model fits (lognormal error) to data of Skalak et al.
Panel (a) shows linearly spaced contours of the objective function
$\ell\left(\mu_{0},\lambda,\diamond\right)$, with a blue diamond showing the
location of the Matlab fmincon parameter fit (with $n^{\diamond}=0.635$), and
a red cross showing the arbitrarily chosen ‘alternative’ parameter choice
(with $n^{\times}=0.634$). Panels (b) and (c) plot the fits of the
Carreau-1972 model to the data of Skalak et al. (1981) for these parameter
tuples shown on linear and logarithmic axes respectively.

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 3, with fitting results in Figures 3b,c.
Fitting over the full range in shear rate (Figure 3) 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 3), 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.

Figure 3: Numerical experimental results: synthetic data generated from, then
fitted to the Carreau-1972 model. (a) Lognormal noise with standard deviation
$0.02$ was used to generate 50 points of synthetic data from the Carreau-1972
flow model with $\mu_{0}=150,\lambda=100,n=0.45$. (b) Cost function, optimized
over $n$ for each $(\mu_{0},\lambda)$ 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).

###  3.2 Cross-1965

As described above, the Cross-1965 model involves the parameters $\mu_{0},\
\lambda,\ n$, and, in addition, an infinite shear rate viscosity
$\mu_{\infty}$. To facilitate comparison with the Carreau-1972 model we will
initially set $\mu_{\infty}=0$, which we will refer to as the ‘Cross-Zero’
model,

$\mu(\dot{\gamma})=\frac{\mu_{0}}{1+(\lambda\dot{\gamma})^{1-n}}.$ (8)

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
$\lambda$ and non-zero shear rate, the constitutive law can be approximated
by,

$\mu(\dot{\gamma})\approx\frac{\mu_{0}}{(\lambda\dot{\gamma})^{1-n}},$ (9)

yielding an infinite family of approximately equivalent parameterizations for
which $\mu_{0}/\lambda^{1-n}$ is constant.

Figure 4: Results of fitting the Skalak et al. data with lognormal error with
our ‘Cross-Zero’ model (equation (8), based on Cross-1965 with
$\mu_{\infty}=0$. Panel (a) shows linearly spaced contours of the objective
function $\ell\left(\mu_{0},\lambda,\diamond\right)$, with a blue diamond
showing the location of the Matlab fmincon parameter fit. Panels (b) and (c)
plot the fits of the Cross-Zero model to the data of Skalak et al. (1981) for
these parameter tuples shown on linear and logarithmic axes respectively.

Having established that this simplified version of the Cross model is also
affected by parameter indeterminancy, we turn our attention to the full four-
parameter 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 5).

Figure 5: Results of fitting the Skalak et al. data with lognormal error with
the Cross-1965 model (equation (1)). Panel (a) shows linearly spaced contours
of the objective function $\ell\left(\mu_{0},\lambda,\diamond\right)$, with a
blue diamond showing the location of the Matlab fmincon parameter fit. Panels
(b) and (c) plot the fits of the Cross-1965 model to the data of Skalak et al.
(1981) for these parameter tuples shown on linear and logarithmic axes
respectively.

###  3.3 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
$\mu_{0}=106,\ \lambda=98.1,\ n=0.635,\ a=7.61$. Nevertheless, there is still
a very elongated ridge in parameter space and associated uncertainty.

Figure 6: Results of fitting the Skalak et al. data with lognormal error with
the Yasuda-1979 model (equation (3)). Panel (a) shows linearly spaced contours
of the objective function $\ell\left(\mu_{0},\lambda,\diamond\right)$, with a
blue diamond showing the location of the Matlab fmincon parameter fit. Panels
(b) and (c) plot the fits of the Cross-Zero model to the data of Skalak et al.
(1981) for these parameter tuples shown on linear and logarithmic axes
respectively.

###  3.4 BCCY-1987

Finally, we consider the most general five-parameter BCCY-1987 model (Figure
7); the best fit parameter tuple (blue diamond) is $\mu_{0}=89.8,\
\mu_{\infty}=3.03,\ \lambda=14.2,\ n=0.339\ a=2.15$ and the cost function
$\ell^{\Diamond}=0.0196$. The fit is excellent, as for the Cross-1965 model,
and as for the Yasuda-1979 model there is a minimum interior to the domain.
However, the indeterminacy is arguably the worst of all models considered,
with a large flat region in the top left corner of the parameter domain
considered. An alternative point chosen at $\mu_{0}=350,\ \lambda=107$ with
optimized $\mu_{\infty}=3.03,\ n=0.334,\ a=0.755$ yields only a marginal
increase in cost function value $\ell^{\times}=0.0214$ (red cross).

Figure 7: Results of fitting the Skalak et al. data with lognormal error with
the BCCY-1987 model (equation (4)). Panel (a) shows linearly spaced contours
of the objective function $\ell\left(\mu_{0},\lambda,\diamond\right)$, with a
blue diamond showing the location of the Matlab fmincon parameter fit, and a
red cross showing the arbitrarily chosen ‘alternative’ parameter choice.
Panels (b) and (c) plot the fits of the BCCY-1987 model to the data of Skalak
et al. (1981) for these parameter tuples shown on linear and logarithmic axes
respectively.

##  4 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
and Cross-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 $\mu_{0}$
and $\lambda$ 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 and BCCY-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 steady-shear
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 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 $\mu_{0}$ and $\lambda$ 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.

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.

## 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

## Acknowledgements

M.T.G. and D.J.S. are supported by the Engineering and Physical Sciences
Research Council award EP/N021096/1. S.D. was supported by the London
Mathematical Society Undergraduate Research Bursary 17-18 13. This work
started with the Multi-scale Biology Study Group, University of Birmingham
(12-15th December 2016), which was jointly funded by POEMS (Predictive
modelling for healthcare technology through maths - EP/L001101/1) and MSB-Net
(UK Multi-Scale Biology Network - BB/M025888/1). The sponsors had no role in
the study design.

## Conflict of interest statement

The authors confirm that there are no financial or personal relationship with
other people or organisations that could inappropriately influence (bias) this
work.

##  Appendix A Indeterminacy of the Carreau-1972 model applied to other
datasets

The parameter fitting to the Carreau-1972 model in Figure 2 is applied to the
combined data sets from Merrill et al. (1963), Cokelet et al. (1963) and Huang
et al. (1973) in Figure 9. Each of these data sets have been extracted from
Ballyk et al. (1994) using GRABIT, Doke (2016).

Figure 9: Likelihood surfaces optimized over $n$ for the Carreau-1972 model
with lognormal error, consistently showing extended flat regions. (a) Merrill
et al., (b) Cokelet et al., (c) Huang et al.

## References

  * Anand & Rajagopal (2017) Anand, M., & Rajagopal, K. R. (2017).  A short review of advances in the modelling of blood rheology and clot formation.  Fluids, 2, 35\.
  * Apostolidis et al. (2015) Apostolidis, A. J., Armstrong, M. J., & Beris, A. N. (2015).  Modeling of human blood rheology in transient shear flows.  Journal of Rheology, 59, 275–298.
  * Ballyk et al. (1994) Ballyk, P., Steinman, D., & Ethier, C. (1994).  Simulation of non-newtonian blood flow in an end-to-side anastomosis.  Biorheology, 31, 565–586.
  * Bird et al. (1987) Bird, R. B., Armstong, R. C., & Hassager, O. (1987).  Dynamics of polymeric liquids Volume 1: fluid mechanics.  (1st ed.).  John Wiley & Sons, Inc.
  * Carreau (1972) Carreau, P. J. (1972).  Rheological Equations from Molecular Network Theories.  Transactions of the Society of Rheology, 16, 99–127.
  * Cokelet et al. (1963) Cokelet, G. R., Merrill, E. W., Gilliland, E. R., Shin, H., Britten, A., & Wells Jr, R. E. (1963).  The rheology of human blood-measurement near and at zero shear rate.  Transactions of the Society of Rheology, 7, 303–317.
  * Cross (1965) Cross, M. M. (1965).  Rheology of non-Newtonian fluids: a new flow equation for pseudoplastic systems.  Journal of Colloid and Interface Science, 20, 417–437.
  * Doke (2016) Doke, J. (2016).  GRABIT.  https://uk.mathworks.com/matlabcentral/fileexchange/7173-grabit.  Accessed: 2018-09-24.
  * Gijsen et al. (1999) Gijsen, F. J. H., van de Vosse, F. N., & Janssen, J. D. (1999).  The influence of the non-Newtonian properties of blood on the flow in large arteries: steady flow in a carotid bifurcation model.  Journal of Biomechanics, 32, 601–608.
  * Huang et al. (1973) Huang, C. R., King, R. G., & Copley, A. L. (1973).  Rheogoniometric studies of whole human blood at shear rates down to 0.0009 sec-1.  Biorheology, 10, 23–28.
  * Kolbasov et al. (2016) Kolbasov, A., Comiskey, P. M., Sahu, R. P., Sinha-Ray, S., Yarin, A. L., Sikarwar, B. S., Kim, S., Jubery, T. Z., & Attinger, D. (2016).  Blood rheology in shear and uniaxial elongation.  Rheologica Acta, 55, 901–908.
  * Leuprecht & Perktold (2001) Leuprecht, A., & Perktold, K. (2001).  Computer simulation of non-newtonian effects on blood flow in large arteries.  Computer Methods in Biomechanics and Biomedical Engineering, 4, 149–163.
  * Merrill et al. (1963) Merrill, E. W., Gilliland, E. R., Cokelet, G., Shin, H., Britten, A., & Wells Jr, R. E. (1963).  Rheology of human blood, near and at zero flow. effects of temperature and haematocrit level.  Biophysical Journal, 3, 199–213.
  * Ostwald (1925) Ostwald, W. (1925).  Ueber die geschwindigkeitsfunktion der viskosität disperser systeme. i.  Kolloid-Zeitschrift, 36, 99–117.
  * Perktold & Rappitsch (1995) Perktold, K., & Rappitsch, G. (1995).  Computer simulation of local blood flow and vessel mechanics in a compliant carotid artery bifurcation model.  Journal of Biomechanics, 28, 845–856.
  * Skalak et al. (1981) Skalak, R., Keller, S. R., & Secomb, T. W. (1981).  Mechanics of blood flow.  Journal of Biomechanical Engineering, 103, 102–115.
  * de Waele (1923) de Waele, A. (1923).  Viscometry and plastometry.  Journal of the Oil & Colour Chemists’ Association, 6, 33–69.
  * Yasuda (1979) Yasuda, K. (1979).  Investigation of the analogies between viscometric and linear viscoelastic properties of polystyrene fluids.  Thesis Massachusetts Institute of Technology.