# On a Finite Field Analogue of van der Corput Differencing

Ikuya Kaneko  The Division of Physics, Mathematics and Astronomy, California
Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA
<EMAIL_ADDRESS>https://sites.google.com/view/ikuyakaneko/

###### Abstract.

We establish nontrivial bounds for general bilinear forms with a given
periodic function, which are thought of as an analogue of van der Corput
differencing for exponential sums. The proof employs Poisson summation,
Cauchy–Schwarz, and the iteration thereof.

###### Key words and phrases:

Van der Corput differencing, Poisson summation, trace functions

###### 2020 Mathematics Subject Classification:

11L40 (primary); 11L20 (secondary)

The author acknowledges the support of the Masason Foundation.

##  1\. Introduction

The notion of exponential sums is one of the most important subjects in
analytic number theory. The problem of determining the optimal size of such
sums is key to understanding many salient problems such as equidistribution of
various elements within domains, counting problems, the distribution of primes
in subsets, and the number of solutions to Diophantine equations, to name a
few. In modern number theory, there are various techniques to obtain their
strong upper bounds; see in particular the pioneering papers of Weyl [Wey16]
and van der Corput [vdC21]. They consider a general exponential sum of the
form

$\sum_{N\leq n<2N}e(f(n)),$

where $f(n)$ is a real-valued function and $e(x)=e^{2\pi ix}$. In order to
prove a nontrivial estimate $N^{1-\delta}$ for some $\delta>0$, van der Corput
begins with replacing the phase $f(x)$ with its derivative $f^{\prime}(x)$ via
the Cauchy–Schwarz inequality, and then applies the Poisson summation formula
so that the dual length is roughly of size $f^{\prime}(N)$. The iteration of
van der Corput differencing leads to a standard $k$-th derivative test, which
now evolves into the theory of exponent pairs. For further details and
historical background, see the book of Graham–Kolesnik [GK91].

This work addresses a finite field analogue of van der Corput differencing for
exponential sums. To enable the ensuing discussion, we introduce some
notation. For a sequence $\alpha=(\alpha_{n})$ supported on a dyadic interval
$n\in[N,2N)$, we write

$\left\lVert\alpha\right\rVert_{2}^{2}\coloneqq\frac{1}{N}\sum_{N\leq
n<2N}|\alpha_{n}|^{2},$

and for a $c$-periodic function $K(m,n;c)$ in $m$ and $n$, we write

$\left\lVert
K\right\rVert_{2,\infty}^{2}\coloneqq\sup_{(m,c)=1}\frac{1}{c}\sum_{x{\@displayfalse\pmod{c}}}|K(m,x;c)|^{2}.$
(1.1)

Note that $\left\lVert\alpha\right\rVert_{2}^{2}$ is normalised by a factor of
$\frac{1}{N}$, which is less standard. Given such a function $K$, we denote

$(\Lambda_{\ell}K)(m,n;c)\coloneqq\frac{1}{\sqrt{c}}\sum_{x{\@displayfalse\pmod{c}}}K(m,x;c)\overline{K(n,x;c)}e\left(-\frac{\ell
x}{c}\right).$

Generally, one expects that if $K$ is a trace function, then so is
$\Lambda_{\ell}K$ and thus $\Lambda_{\ell}K\ll 1$. We then define a bilinear
form by

$\mathcal{S}_{K}(M,N;c)\coloneqq\sum_{M\leq m<2M}\sum_{N\leq
n<2N}\alpha_{m}\beta_{n}K(m,n;c).$ (1.2)

The main result of this paper is the following.

######  Theorem 1.1.

Let $K(m,n;c)$ be a $c$-periodic function in $m$ and $n$. Let
$\alpha=(\alpha_{m})$ and $\beta=(\beta_{n})$ traverse sequences supported on
$m\in[M,2M)$ and $n\in[N,2N)$, respectively. Suppose that $M\asymp N$. Then we
have that

$\mathcal{S}_{K}(M,N;c)\ll\sqrt{cMN}\left\lVert\alpha\right\rVert_{2}\left\lVert\beta\right\rVert_{2}\sup_{|\ell_{1}|,\ldots,|\ell_{k}|\ll\frac{c}{N}}\bigg{(}\left\lVert
K\right\rVert_{2,\infty}+\left\lVert\Lambda_{\ell_{1}}K\right\rVert_{2,\infty}^{\frac{1}{2}}+\cdots\\\
+\left\lVert\Lambda_{\ell_{k-1}}\cdots\Lambda_{\ell_{1}}K\right\rVert_{2,\infty}^{\frac{1}{2^{k-1}}}+M^{\frac{1}{2^{k+1}}}\left\lVert\alpha\right\rVert_{2}^{-\frac{1}{2^{k}}}\left\lVert\Lambda_{\ell_{k}}\cdots\Lambda_{\ell_{1}}K\right\rVert_{2,\infty}^{\frac{1}{2^{k}}}\bigg{)}.$
(1.3)

The proof of Theorem 1.1 features Poisson summation, Cauchy–Schwarz, and the
iteration ad infinitum thereof. Despite the softness of our input, it has a
square-root improvement over the trivial bound of size $MN$ (if we assume that
all the norms are bounded). Furthermore, Theorem 1.1 keeps the coefficients
$(\alpha_{m})$ and $\beta_{n}$ as general as possible, which would facilitate
potential arithmetic implications for the subconvexity problem for families of
$L$-functions for instance. In particular, a version of the bilinear form
(1.2) often arises when one utilises the delta method to separate
oscillations, in the terminology of Munshi [Mun18]. In this respect, an
interesting case is when $K$ is a hyper-Kloosterman sum, in which case (1.2)
is studied by Kowalski–Michel–Sawin [KMS17] and Kerr et al. [KSWX23]. Note
that these iterative ideas are present in [ALM22, HB78, KL23], but not written
down in a general format.

######  Remark 1.2.

The assumption $M\asymp N$ in Theorem 1.1 is imposed to decrease the
complexity of the argument, and in fact our strategy should be capable of
dealing with a wider regime.

### Acknowledgements

The author thanks Maksym Radziwiłł for suggesting underlying ideas in this
paper, and Wing Hong Leung for productive discussions.

##  2\. Preliminaries

For $n\in\mathbb{N}$ and an integrable function
$w\colon\mathbb{R}^{n}\to\mathbb{C}$, we denote its Fourier transform by

$\widehat{w}(\xi)\coloneqq\int_{\mathbb{R}^{n}}w(t)e(-\langle
t,\xi\rangle)dt,$

where $\langle\cdot,\cdot\rangle$ denotes the standard inner product on
$\mathbb{R}^{n}$. Furthermore, if $c\in\mathbb{N}$ and
$K\colon\mathbb{Z}\to\mathbb{C}$ is a periodic function of period $c$, then
its discrete Fourier transform $\widehat{K}$ is again the periodic function of
period $c$ given by

$\widehat{K}(n)\coloneqq\sum_{a{\@displayfalse\pmod{c}}}K(a)e\left(-\frac{an}{c}\right).$

We invoke a version of the Poisson summation formula with a $c$-periodic
function involved.

######  Lemma 2.1 (Fouvry–Kowalski–Michel [FKM15, Lemma 2.1]).

For any $c\in\mathbb{N}$, any $c$-periodic function $K$, and any even smooth
function $V$ compactly supported on $\mathbb{R}$, we have that

$\sum_{n=1}^{\infty}K(n)V(n)=\frac{1}{c}\sum_{n\in\mathbb{Z}}\widehat{K}(n)\widehat{V}\left(\frac{n}{c}\right).$

##  3\. Proof of Theorem 1.1

To eschew excessive complications, we can assume without loss of generality
that $\alpha_{m}$ is such that $\alpha_{m}\neq 0$ implies $m$ is prime and
$(m,c)=1$. We further assume that $|K(m,n;c)|\leq 1$ for $(m,n)=1$, and
$|K(m,n;c)|\leq\sqrt{c}$ otherwise. Throughout this section, the notation
$m\sim M$ serves as a shorthand for $M\leq m<2M$.

First of all, the Cauchy–Schwarz inequality implies that there exists a
compactly supported smooth function $V\in C_{c}^{\infty}([1,2])$ such that

$|\mathcal{S}_{K}(M,N;c)|\leq
N\left\lVert\beta\right\rVert_{2}\left(\frac{1}{N}\sum_{n}\left|\sum_{m\sim
M}\alpha_{m}K(m,n;c)\right|^{2}V\left(\frac{n}{N}\right)\right)^{\frac{1}{2}}.$

Squaring out and applying Poisson summation (Lemma 2.1) to the sum over $n$,
we obtain

$c^{\frac{1}{4}}\sqrt{N}\left\lVert\beta\right\rVert_{2}\left(\frac{N}{c}\sum_{\ell}\sum_{m_{1},m_{2}\sim
M}\alpha_{m_{1}}\alpha_{m_{2}}(\Lambda_{\ell}K)(m_{1},m_{2};c)\widehat{V}\left(\frac{\ell}{c/N}\right)\right)^{\frac{1}{2}}.$

Repeated integration by parts ensures an arbitrary saving unless
$|\ell|\ll\frac{c}{N}$, thus we arrive at

$c^{\frac{1}{4}}\sqrt{N}\left\lVert\beta\right\rVert_{2}\left(\frac{N}{c}\sum_{|\ell|\ll\frac{c}{N}}\sum_{m_{1},m_{2}\sim
M}\alpha_{m_{1}}\alpha_{m_{2}}(\Lambda_{\ell}K)(m_{1},m_{2};c)\right)^{\frac{1}{2}}.$

We decompose this expression according as $m_{1}=m_{2}$ or $m_{1}\neq m_{2}$.
Upon bounding the sum over $\ell$ trivially, the contribution of the diagonal
terms $m_{1}=m_{2}$ boils down to

$\sqrt{cMN}\left\lVert\beta\right\rVert_{2}\left(\frac{1}{M}\sum_{m\sim
M}|\alpha_{m}|^{2}\cdot\frac{1}{c}\sum_{x{\@displayfalse\pmod{c}}}|K(m,x;c)|^{2}\right)^{\frac{1}{2}}.$

Using the notation (1.1), we may estimate the diagonal contribution as

$\sqrt{cMN}\left\lVert\alpha\right\rVert_{2}\left\lVert\beta\right\rVert_{2}\left\lVert
K\right\rVert_{2,\infty}.$

When $m_{1}\neq m_{2}$, the condition $(m_{1},m_{2})=1$ holds automatically.
Hence we may estimate the off-diagonal contribution as

$c^{\frac{1}{4}}\sqrt{N}\left\lVert\beta\right\rVert_{2}\mathcal{S}_{\Lambda_{\ell}K}(M,M,c)^{\frac{1}{2}}$

for some $|\ell|\ll\frac{c}{N}$. To summarise, there holds

$\mathcal{S}_{K}(M,N;c)\ll\sqrt{cMN}\left\lVert\alpha\right\rVert_{2}\left\lVert\beta\right\rVert_{2}\left\lVert
K\right\rVert_{2,\infty}+c^{\frac{1}{4}}\sqrt{N}\left\lVert\beta\right\rVert_{2}\mathcal{S}_{\Lambda_{\ell}K}(M,M,c)^{\frac{1}{2}}.$

We now start iterating the above process deducing

$\mathcal{S}_{K}(M,N;c)\ll\sqrt{cMN}\left\lVert\alpha\right\rVert_{2}\left\lVert\beta\right\rVert_{2}\left\lVert
K\right\rVert_{2,\infty}+\sqrt{cMN}\left\lVert\alpha\right\rVert_{2}\left\lVert\beta\right\rVert_{2}\left\lVert\Lambda_{\ell_{1}}K\right\rVert_{2,\infty}^{\frac{1}{2}}\\\
+c^{\frac{3}{8}}\sqrt{N}M^{\frac{1}{4}}\left\lVert\alpha\right\rVert_{2}^{\frac{1}{2}}\left\lVert\beta\right\rVert_{2}\mathcal{S}_{\Lambda_{\ell_{2}}\Lambda_{\ell_{1}}K}(M,M,c)^{\frac{1}{4}}.$

Iterating $k-1$ times for sufficiently large $k\in\mathbb{N}$ implies

$\mathcal{S}_{K}(M,N;c)\ll\sqrt{cMN}\left\lVert\alpha\right\rVert_{2}\left\lVert\beta\right\rVert_{2}\left(\left\lVert
K\right\rVert_{2,\infty}+\left\lVert\Lambda_{\ell_{1}}K\right\rVert_{2,\infty}^{\frac{1}{2}}+\cdots+\left\lVert\Lambda_{\ell_{k-1}}\cdots\Lambda_{\ell_{1}}K\right\rVert_{2,\infty}^{\frac{1}{2^{k-1}}}\right)\\\
+c^{\frac{1}{2}-\frac{1}{2^{k+1}}}\sqrt{N}M^{\frac{1}{2}-\frac{1}{2^{k}}}\left\lVert\alpha\right\rVert_{2}^{1-\frac{1}{2^{k-1}}}\left\lVert\beta\right\rVert_{2}\mathcal{S}_{\Lambda_{\ell_{k}}\cdots\Lambda_{\ell_{1}}K}(M,M,c)^{\frac{1}{2^{k}}}.$
(3.1)

To conclude the iteration, we employ the trivial bound

$|\mathcal{S}_{K}(M,N;c)|\leq\sqrt{MN}\left\lVert\alpha\right\rVert_{2}\left\lVert\beta\right\rVert_{2}\left(\mathop{\sum\sum}_{\begin{subarray}{c}m\sim
M,\,n\sim N\\\
(m,n)=(mn,c)=1\end{subarray}}|K(m,n;c)|^{2}\right)^{\frac{1}{2}}.$

Applying Poisson summation (Lemma 2.1) in either of the variables shows

$\sum_{n\sim
N}|K(m,n;c)|^{2}\ll\frac{N}{c}\sum_{|\ell|\ll\frac{c}{N}}\sum_{x{\@displayfalse\pmod{c}}}|K(m,x;c)|^{2}e\left(-\frac{\ell
x}{c}\right)\ll c\left\lVert K\right\rVert_{2,\infty}^{2}.$

Hence

$\mathcal{S}_{K}(M,N;c)\ll\min(\sqrt{M},\sqrt{N})\sqrt{cMN}\left\lVert\alpha\right\rVert_{2}\left\lVert\beta\right\rVert_{2}\left\lVert
K\right\rVert_{2,\infty}.$

Inserting this estimate into (3.1) completes the proof of Theorem 1.1.

## References

  * [ALM22] Keshav Aggarwal, Wing Hong Leung, and Ritabrata Munshi, _Short Second Moment Bound and Subconvexity for $\mathrm{GL}(3)$ $L$-Functions_, arXiv e-prints (2022), 58 pages.
  * [FKM15] Étienne Fouvry, Emmanuel Kowalski, and Philippe Gabriel Michel, _On the Exponent of Distribution of the Ternary Divisor Function_ , Mathematika 61 (2015), no. 1, 121–144. MR 3333965
  * [GK91] Sidney West Graham and Grigori Kolesnik, _Van der Corput’s Method of Exponential Sums_ , London Mathematical Society Lecture Note Series, vol. 126, Cambridge University Press, Cambridge, 1991. MR 1145488
  * [HB78] Roger Heath-Brown, _The Twelfth Power Moment of the Riemann-Function_ , The Quarterly Journal of Mathematics 29 (1978), no. 116, 443–462. MR 517737
  * [KL23] Ikuya Kaneko and Wing Hong Leung, _Non-Linear Additive Twists of $\mathrm{GL}_{3}$ Hecke Eigenvalues_, preprint (2023), 26 pages.
  * [KMS17] Emmanuel Kowalski, Philippe Gabriel Michel, and William French Sawin, _Bilinear Forms with Kloosterman Sums and Applications_ , Annals of Mathematics 186 (2017), no. 2, 413–500. MR 3702671
  * [KSWX23] Bryce Kerr, Igor E. Shparlinski, Xiaosheng Wu, and Ping Xi, _Bounds on Bilinear Forms with Kloosterman Sums_ , Journal of the London Mathematical Society 108 (2023), no. 2, 578–621. MR 4626720
  * [Mun18] Ritabrata Munshi, _The Subconvexity Problem for $L$-Functions_, Proceedings of the International Congress of Mathematicians—Rio de Janeiro 2018. Volume II. Invited Lectures. ICM 2018, August 1–9, 2018 (Boyan Sirakov, Paulo Ney de Souza, and Marcelo Viana, eds.), World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2018, pp. 363–376. MR 3966770
  * [vdC21] Johannes Gaultherus van der Corput, _Zahlentheoretische Abschätzungen_ , Mathematische Annalen 84 (1921), no. 1-2, 53–79. MR 1512020
  * [Wey16] Hermann Klaus Hugo Weyl, _Über die Gleichverteilung von Zahlen mod. Eins_ , Mathematische Annalen 77 (1916), no. 3, 313–352. MR 1511862


# Effects of a Mixed Reality Headset on the Delay of Visually Evoked
Potentials

Víctor Manuel Hidalgo  Research and Development Department, DTS SpA, National
Aeronautical Company of Chile  Corresponding author: Victor Manuel Hidalgo,
<EMAIL_ADDRESS>Carlos Andrés Bazaes  Research and Development Department,
DTS SpA, National Aeronautical Company of Chile  Department of Biology, School
of Sciences, University of Chile  Juan-Carlos Letelier  Department of Biology,
School of Sciences, University of Chile

###### Abstract

Virtual and mixed reality (VR, MR) technologies offer a powerful solution for
on-the-ground flight training curricula. While these technologies offer safer
and cheaper instructional programs, it is still unclear how they impact
neuronal brain dynamics. Indeed, MR simulations engage students in a strange
mix of incongruous visual, somatosensory and vestibular sensory input.
Characterizing brain dynamics during MR simulation is important for
understanding cognitive processes during virtual flight training. To this end,
we studied the delays introduced in the neuronal stream from the retina to the
visual cortex when presented with visual stimuli using a Varjo-XR3 headset. We
recorded cortical visual evoked potentials (VEPs) from 6 subjects under two
conditions. First, we recorded normal VEPs triggered by short flashes. Second,
we recorded VEPs triggered by an internal image of the flashes produced by the
Varjo-XR3 headset. All subjects had used the headset before and were familiar
with immersive experiences. Our results show mixed-reality stimulation imposes
a small, but consistent, 4 [ms] processing delay in the N2-VEP component
during MR stimulation as compared to direct stimulation. Also we found that
VEP amplitudes during MR stimulation were also decreased. These results
suggest that visual cognition during mixed-reality training is delayed, not
only by the unavoidable hardware/software processing delays of the headset and
the attached computer, but also by an extra biological delay induced by the
headset’s limited visual display in terms of image intensity and contrast. As
flight training is a demanding task, this study measures visual signal latency
to better understand how MR affects the sensation of immersion.

##  1 Introduction

Mixed reality simulators are a promising solution for on-the-ground flight
training given the feature that allows interaction with both virtual and real
environments. Applying mixed reality is still relatively new, but using this
approach for training has great potential in multiple areas and has the added
value of reduced risk and cost. However, the effects of virtual reality
experiences on the brain are not fully understood. There is some research on
the use of EEG in Virtual Reality [6], but how mixed reality alters the visual
experience of subjects – thus affecting training – remains an open question.
Measuring visual evoked potentials (VEPs) is one approach that has been
heavily researched to recognize the effects of diverse variables such as
stimulus intensity [3], anomalies in the subject’s eyes [7], and stimulus
adaption [5]. This type of measurement was thus chosen to investigate the
effects of a mixed-reality Head Mounted Display (HMD) on visual processing. We
found that visual evoked potentials are significantly delayed when elicited by
a mixed reality HMD, suggesting that there is an extra biological delay in
experience when using mixed reality systems.

Figure 1: Mixed reality flight simulator MUPUN-X. Mixed reality has the
capacity to show both the real environment and a virtual scenario at the same
time. In this example, the pilot can see the cockpit while a flight simulation
is displayed through the headset on the green screen.

##  2 Methods

###  2.1 Subjects

For this study, 6 subjects (5 men, 1 woman, ages=(21,22,23,23,26,64)) were
recruited from the University of Chile. The experiment was approved by the
”Comité de Ética de la Facultad de Ciencias de la Universidad de Chile”
(Permit ID: 2225-FCS-UCH), and all subjects signed a formal consent form
before any experiment was conducted.

###  2.2 Experiments

Each subject was seated at 15 [cm] in front of a light flash stimulator made
by a row of 5 bright white LEDs. An Arduino Nano controlled this light source
and was set to generate 100 [ms] stimuli separated by an interval uniformly
distributed between 3 and 7 [s]. Visual Evoked Potentials (VEPs) were recorded
using a single-channel EEG from F0-O2 under 2 conditions. In condition 1
(_c1_), subjects’ VEPs were recorded when they were directly stimulated
without using the HMD. In condition 2 (_c2_) subjects used the mixed reality
headset Varjo XR-3 with the resolution set to ”very-high” (1984x1984
1976x1696, 65 pixels for degree) and a 90Hz refresh rate. The hardware
controlling this headset consisted of a Notebook ORIGIN EON17-X with a
I9-11900K processor, a NVIDIA GeForce RTX 3080 16GB GDDR6 dedicated graphic
card, and 64 GB of ram memory, under the windows 11 Pro operating system.

To record the stimulus timing we used a light sensor. In condition _c1_ (no
HMD) the sensor directly faced the light stimulator. In condition _c2_ (i.e.
use of the Varjo HMD) the light sensor was located inside the HMD to measure
the time of appearance of the flash images in the inner uOLED Varjo displays.
Subjects were instructed to look at the light stimulator for a 200 [s] time
period, take a short rest, and then repeat for a new period until 160 VEPs
were collected.

Figure 2: VEP recording setup. A LED light source is placed 15 [cm] away from
the subject and visually evoked potentials are measured in 2 conditions. Left
panel (_c1_ ): VEPs are recorded using a direct visual stimulus. Right panel
_c2_ : VEPs are recorded while subjects use a head mounted mixed reality
display; in this case, VEPs are elicited by the flash images on the HMD uOLED
displays.

###  2.3 Analysis

All analyses were conducted with IGOR 9 PRO. For each VEP recording, we
isolated 100 [ms] before and 300 [ms] after the stimulus from the EEG
recording. Those recordings with ocular movement or with very distorted
signals were set aside. From the remaining data, 75 random VEPs were selected
for each subject and condition. These recordings were manually analyzed to
identify the N2 peak at around 75 [ms] [4, 1]. To obtain valid statistics to
compare the delays of N2 peaks, we applied a resampling test (n=20) for 1000
instances for conditions and subjects [2].

##  3 Results

The distribution of all visual evoked potentials (VEPs) for each condition and
each subject is shown in Figure 3. In general, the distribution of condition 2
has a lower amplitude and a greater standard deviation. Furthermore delays
introduced by the normal processing of the visual image (i.e. initial capture
of the external image plus image processing in the host computer and the final
download onto internal displays) are around 35 [ms] (data not shown).

Figure 3: Summary of Visual Evoked Potentials for each subject and each
condition. Each panel represents the collection of VEPs recorded for each
subject (red=condition _c1_ , grey=condition _c2_ , heavy traces represent the
ensemble average, and the diffuse colored cloud represents the standard
deviation for a given time). The dotted horizontal line indicates the 0 value
for the corresponding VEPs for conditions _c1_ and _c2_. The dashed line
represents stimulus onset time (_c1_ : LEDs ON, _c2_ : internal image ON).
Arrows indicate the N2 peak for each subject

To visualize if there is a difference in the perception of the stimuli during
use of the HMD, we compared the delays of the N2 peak of each VEP condition
_c1_ and _c2_ , then applied a resampling test for each condition and subject
(data is summarized in Figure 4). We then applied a Two-Sample t-Test with a
Bonferroni correction $\alpha=0.05$ for each subject to compare the latency
data under the two conditions and reject the null hypothesis ($H_{0}:$
Latencies are the same).

Figure 4: Resampling distributions of N2 latency for each subject and
experimental condition. Each panel shows the histogram for the mean of 20
resampled N2 latency samples (1000 resampling instances) for each condition
(Condition 1 is represented in red and Condition 2 in black). A two-tailed
two-sample t-test at $\alpha=0.05$ with Bonferroni correction rejects the null
hypothesis $\mu_{C1}=\mu_{C2}$ for all subjects.

##  4 Discussion

In this work, we obtained EEG recordings to explore the effects of a mixed-
reality HMD on the delay of visual evoked potentials (VEPs) triggered by short
flashes. We showed that a high performance HMD lens introduces a small delay
in visual processing. We found that the N2 peak latency was significantly
delayed, by around 4 ms when the HMD was used. This difference in N2 could be
related with a reduction in the intensity or in the contrast of the HMD images
of the external flashes [3]. Mixed reality is a promising technology for
virtual flight training but its effects on visual and especially cortical
processing are not fully understood. We hypothesized that the small change in
latency introduced by HMD is related to the drastic change occurring with
respect to the intensity and contrast of the real image vis-a-vis the image
displayed. As VEPs reflect the activity of the underlying neuronal
populations, it is not surprising that VEPs triggered by stimuli presented in
the VARJO uOLED internal displays have less amplitude and exhibit a slightly
larger latency [3]. This result suggests that mixed reality training imposes a
small, but consistent, delay in activating the visual cortex. In addition,
each HMD introduces a digital processing delay due to the unavoidable steps of
image capture; in the case of VARJO XR-3 this delay is around 35 [ms].
Overall, our result shows that a high performance HMD introduces delays due to
hardware/software issues and neuronal dynamics. The neuronal delay albeit not
negligible is small enough to enable the operational closure needed to produce
the sensation of immersion. If we measure the observed delay in terms of frame
rate (90Hz is equivalent to a 11.1 [ms]), it would induce a maximum delay
equivalent to a single video frame. In any case, sophisticated augmented
reality projects demand that these optic, visual and neuronal parameters are
measured and understood in order to produce a vivid, useful immersive learning
system.

## References

  * [1] American Clinical Neurophysiology Society  “Guideline 9B: Guidelines on visual evoked potentials”  In _J Clin Neurophysiol_ 23.2, 2006, pp. 138–156
  * [2] Anthony Davison and D. Hinkley  “Bootstrap Methods and Their Application”  In _Journal of the American Statistical Association_ 94, 1997  DOI: 10.2307/1271471
  * [3] R.G. DeVoe, H. Ripps and H.G. Vaughan  “Cortical responses to stimulation of the human fovea”  In _Vision Research_ 8.2, 1968, pp. 135–147  DOI: https://doi.org/10.1016/0042-6989(68)90002-3
  * [4] J Vernon Odom et al.  “Visual evoked potentials standard (2004)”  In _Doc Ophthalmol_ 108.2, 2004, pp. 115–123
  * [5] W.H. Ridder and S. Nusinowitz  “The visual evoked potential in the mouse—Origins and response characteristics”  In _Vision Research_ 46.6, 2006, pp. 902–913  DOI: https://doi.org/10.1016/j.visres.2005.09.006
  * [6] Jan-Philipp Tauscher et al.  “Immersive EEG: Evaluating Electroencephalography in Virtual Reality”  In _2019 IEEE Conference on Virtual Reality and 3D User Interfaces (VR)_ , 2019, pp. 1794–1800  DOI: 10.1109/VR.2019.8797858
  * [7] Xingyu Yang et al.  “F105. The influence of myopia on visual evoked potential” Abstracts of the 31st International Congress of Clinical Neurophysiology (ICCN) of the IFCN, May 1–6, 2018, Washington, DC, USA  In _Clinical Neurophysiology_ 129, 2018, pp. e105–e106  DOI: https://doi.org/10.1016/j.clinph.2018.04.268


# Systematic evaluation of the population-level effects of alternative
treatment strategies on the basic reproduction number

Dmitry Gromov<EMAIL_ADDRESS>Ingo Bulla<EMAIL_ADDRESS>Ethan O.
Romero-Severson<EMAIL_ADDRESS>Faculty of Applied Mathematics and Control
Processes, Saint Petersburg State University,
St. Petersburg, Russia  Institut für Mathematik und Informatik, Universität
Greifswald, Walther-Rathenau-Straße 47, 17487 Greifswald, Germany  Université
Perpignan Via Domitia, IHPE UMR 5244, CNRS, F-66860 Perpignan, France
Theoretical Biology and Biophysics Group, Los Alamos National Laboratory,
Los Alamos, New Mexico, USA

###### Abstract

An approach to estimate the influence of the treatment-type controls on the
basic reproduction number, $R_{0}$, is proposed and elaborated. The presented
approach allows one to estimate the effect of a given treatment strategy or to
compare a number of different treatment strategies on the basic reproduction
number. All our results are valid for sufficiently small values of the
control. However, in many cases it is possible to extend this analysis to
larger values of the control as was illustrated by examples.

###### keywords:

Disease control, Basic reproduction number , Treatment , Compartmental
epidemic model , Next generation matrix

##  1 Introduction

The basic reproduction number ($R_{0}$), i.e. the number of infections
generated by an infected person in a fully susceptible population, is a key
determinant of the dynamics of an infectious disease. Interpretation of
$R_{0}$ is complex and is covered in a number of excellent papers, see, e.g.,
[28, 15, 17] as well as the books [11, 10]. For the purposes of this paper, we
note that $R_{0}$ can be interpreted as a threshold parameter: if $R_{0}<1$
then the disease can not sustain itself in the population and will eventually
die out. That is, if our goal is to push a disease towards extinction, then we
can categorize a population-level intervention as ‘effective’ if it reduces
$R_{0}$ from where it was before the intervention was enacted.

#### Computation of $R_{0}$

$R_{0}$ can be computed or estimated from partially empirical considerations,
for instance the “survival function”-based approach [16] that assumes
knowledge of both the survival probability $F(t)$ and the infectivity $b(t)$
of an individual as functions of time. In addition to being difficult to
implement, such approaches can hardly be used for computing $R_{0}$ in the
case where there are more than one group of infected. This issue was addressed
in [12], where an elegant approach to computing $R_{0}$ was proposed, termed
the next generation matrix (NGM) method. This method boils down to computing
$R_{0}$ as the spectral radius of a specially constructed matrix. Later, this
method was detailed and substantiated in a number of papers, see, e.g., [28]
and [8]. The power of the NGM approach lies in its universality. It can be
applied to any population balance model as long as it satisfies a number of
very natural assumptions. The result is based on certain properties of
positive and inverse positive matrices, in particular M-matrices [3, 19].

#### Contribution

We aim at developing this approach in that we extend the notion of the basic
reproduction number to a class of controlled disease propagation models. Much
of applied theoretical epidemiology focuses on how treatment policies,
intervention designs, and novel treatments will impact the burden of a given
disease by considering alternative scenarios. There have been a large number
of papers aimed at evaluating the efficacy of different treatment schemes
using various methods, from numerical simulation to optimal control. We
mention [7, 6, 21, 4, 2, 26, 23, 27, 14] for a short list of related research.
Modeling efforts have even attempted to measure the effects of control
programs in historical epidemics [5]. However, most modeling efforts do not
attempt to measure intervention effects systematically, i.e. by considering
the joint effects of model parameters on the relative efficacy of alternative
interventions. In this paper we address treatment programs, i.e., the
intervention strategies that result in moving infected individuals either into
a different group of infected (with different biological or behavioral
characteristics) or into a group of not-infected, which can correspond to
susceptible, recovered or any other group that consists of not contagious
individuals.

When applied to one or several groups of infected individuals, the action of a
treatment strategy can be described by a parameter $u$ that corresponds to a
fraction of potentially eligible individuals that are administered the
treatment during a unit of time. Typically, this fraction is rather small. Our
first result consists in defining the notion of a controlled reproduction
number $R_{0}(u)$ which explicitly contains $u$ as a parameter. The next
result allows to estimate the action of the particular treatment strategy on
the value of the controlled reproduction number $R_{0}(u)$. We formulate
conditions under which the application of a treatment strategy leads to the
reduction in $R_{0}(u)$. This allows for evaluating the efficacy of a devised
treatment strategy as well as for comparing the efficacy of different
strategies. The obtained results are illustrated by a number of examples.

The paper is organized as follows: Section 2 describes a general compartmental
epidemic model and briefly introduces the NGM method for computing the basic
reproduction number $R_{0}$. The main result of the paper along with a
discussion on further extensions and ramifications of the developed approach
are presented in Sec. 3. The obtained results are illustrated by a number of
examples in Sec. 4. The paper is concluded with a discussion section. Finally,
there are two appendices containing necessary technical information.

##  2 Epidemics dynamics

###  2.1 A disease propagation model

When modeling the process of disease propagation in a heterogeneous
population, the standard procedure consists in dividing the total population
into a number of classes (compartments) according to some criteria relevant to
the disease transmission: disease status, ability to contract disease,
behavioral (contact) pattern and so on. All individuals within a compartment
are assumed to be identical in their evolution. Thus we can consider only the
number of individuals within each compartment, denoted by $x_{i}$, where
$i=1,\ldots,n$ is the number of the compartment.

#### General formulation

The dynamics of these groups include transitions between groups and the in-
and outflows associated with these groups. In the following it is assumed that
the total inflow is constant over time and does not depend on the population
size while the outflow depends on both the size and the structure of the
population. The evolution of the $i$th state is thus described by the
following compartmental DE:

$\dot{x}_{i}=\Phi_{i}(x)=w_{i}+\sum\limits_{i\neq
j}\big{(}a_{ij}(x)-a_{ji}(x)\big{)}-a_{i}x_{i},$ (1)

where $x\in\mathbb{R}^{n}_{\geq 0}$ is the state, $w_{i}$ is the constant
inflow, $a_{i}x_{i}$ is the outflow from the $i$-th compartment, and
$a_{ij}(x)$, $i\neq j$, is the flow rate from the $j$-th to the $i$-th
compartment.

The flow rate functions are assumed to be $C^{\infty}$ for all
$x\in\mathbb{R}^{n}_{\geq 0}$. Furthermore, all flow rates have to satisfy the
following properties:

$\displaystyle w_{i}\geq 0,a_{i}\geq 0,a_{ji}(x)\geq 0,$ $\displaystyle\forall
x\in\mathbb{R}^{n}_{\geq 0},i,j=1,\dots,n.$ (2a) $\displaystyle
x_{i}=0\implies a_{ji}(x)=0$ $\displaystyle\forall i,j=1,\dots,n.$ (2b)

The property (2b) implies that there is no outflow from an empty compartment.

#### Compartmental epidemic model

All state variables are divided into two groups depending on whether the
respective compartment corresponds to the infectious individuals (that is
those able to transmit infection) or to the non-infectious individuals
(regardless of whether they were infected earlier or not).

Let there be $l$ state variables representing the “infectious” compartments.
We rearrange the states in the way that the vector of the state variables
takes the form $x^{\top}=[x^{\top}_{I},x_{N}^{\top}]$, where
$x^{\top}_{I}=[x_{1},\dots,x_{l}]$ and $x^{\top}_{N}=[x_{l+1},\dots,x_{n}]$
are the states associated with “infectious” and “non-infectious” compartments.
We will refer to the respective compartments as the $I$\- and
$N$-compartments. Further, we write $\Phi(x)$ as the sum of two vector-valued
functions

$\Phi(x)=\Phi^{I}(x)+\left[w+\Phi^{in}(x)-\Phi^{out}(x)\right]=\Phi^{I}(x)+\Phi^{C}(x),$
(3)

where $\Phi_{i}^{out}(x)=\sum\limits_{j\neq i}a_{ji}(x)+a_{i}x_{i}$, and

$\Phi^{I}(x)=\begin{bmatrix}\bigg{(}\sum\limits_{j>l}a_{ij}(x)\bigg{)}_{i=1,\dots,l}\\\\[20.0pt]
{\mathbf{0}}\end{bmatrix},\enskip\Phi^{in}(x)=\begin{bmatrix}\bigg{(}\sum\limits_{\begin{subarray}{c}j\leq
l\\\ j\neq i\end{subarray}}a_{ij}(x)\bigg{)}_{i=1,\dots,l}\quad\\\\[10.0pt]
\bigg{(}\sum\limits_{j\neq i}a_{ij}(x)\bigg{)}_{i=l+1,\dots,n}\end{bmatrix}.$

Here, $\Phi_{i}^{I}$ is the rate at which new infections occur in the $i$-th
compartment (new infection means that the flow is to be from an
$N$-compartment), $\Phi^{in}(x)$ is the vector of inter-compartmental inflow
rates not related to new infections, and $\Phi^{out}(x)$ are the outflow rates
from the respective compartments. Thus, $\Phi^{C}(x)$ represents the totality
of all flows not related to new infections. Note that $\Phi^{I}_{i}(x)=0$ for
all $i>l$ as there is no inflow of infected into the $N$-compartments.

The components of the respective vectors $\Phi^{(\cdot)}(x)$ satisfy the
following properties (note that (4a) and (4b) can be derived from (2)):

$\displaystyle\Phi_{i}^{I}(x),\Phi_{i}^{in}(x),\Phi_{i}^{out}(x)\geq 0,$
$\displaystyle\forall i=1,\dots,n$ (4a)
$\displaystyle\\{x_{i}=0\\}\implies\\{\Phi_{i}^{out}(x)=0\\},$
$\displaystyle\forall i=1,\dots,n$ (4b)
$\displaystyle\\{x_{I}=\mathbf{0}\\}\implies\\{\Phi_{i}^{I}(x)=0,\Phi_{i}^{in}(x)=0\\},$
$\displaystyle\forall i=1,\dots,l.$ (4c)

The latter property implies that no new infections occur in a totally healthy
population.

#### Assumptions

We conclude this section by making a number of epidemiologically motivated
assumptions.

  * A1.

There is no steady inflow into the $I$-compartments, i.e., $w_{i}=0$,
$i=1,\dots,l$.

  * A2.

There exists a positive constant $\mu$ such that $a_{i}\geq\mu>0$ for all
$i=1,\dots,n$. We refer to $\mu$ as the baseline mortality rate.

###  2.2 Stability of a disease-free equilibrium

Let $\mathcal{X}_{DF}\subset\mathbb{R}^{n}_{\geq 0}$ be the set of disease-
free states, i.e., $\mathcal{X}_{DF}=\\{x\in\mathbb{R}^{n}_{\geq
0}|x_{i}=0,i=1,\dots,l\\}$.

######  Definition 2.1.

Let $x^{*}$ be an equilibrium state, i.e., $\Phi(x^{*})=0$. Then $x^{*}$ is
said to be a disease-free equilibrium (DFE) if $x^{*}\in\mathcal{X}_{DF}$.
Otherwise, $x^{*}$ is referred to as an endemic infection equilibrium.

The system (1) is locally stable at a DFE $x^{*}$ if the linearized model
$\dot{\tilde{x}}=A\tilde{x}$ is asymptotically stable. The latter holds if the
spectrum of the structure matrix $A$ contains only the eigenvalues with
negative real part, [9]. Such matrices are said to be Hurwitz.

######  Lemma 2.1 ([28]).

Let $x^{*}$ be a DFE. The Jacobian matrices
$A^{I}=D\Phi^{I}(x)\big{|}_{x=x^{*}}$ and
$A^{C}=D\Phi^{C}(x)\big{|}_{x=x^{*}}$ have the following form:

$A^{I}=\begin{bmatrix}A^{I}_{11}&0\\\\[10.0pt] 0&0\end{bmatrix},\quad
A^{C}=\begin{bmatrix}A^{C}_{11}&0\\\\[10.0pt]
A^{C}_{21}&A^{C}_{22}\end{bmatrix},$

where $A^{I}_{11}$ is an $[l\times l]$ matrix with non-negative elements,
$A^{I}_{11}\succeq 0$ and $A^{C}_{11}$ is an $[l\times l]$ matrix with non-
positive diagonal elements and non-negative off-diagonal elements111See
Appendix B for the explanation of the used notation..

At this point, we make one more assumption regarding the behavior of the
system in the absence of the infection.

  * A3.

$A^{C}$ is Hurwitz.

This assumption effectively implies that the DFE is asymptotically stable
provided the virus or whatever source of the infection has lost its
contagiousness hence, no new infections occur.

Since the block matrix $A^{C}_{22}$ is stable at a DFE, stability of the whole
system is determined by the sum $A^{I}_{11}+A^{C}_{11}$. Following the common
convention we write $F=A^{I}_{11}$ and $V=-A^{C}_{11}$. Hence, the stability
of a DFE is determined by the stability of the matrix $F-V$.

######  Theorem 2.2 ([28]).

The matrix $F-V$ is Hurwitz if and only if $\rho(FV^{-1})<1$.

We conclude this section with the following two definitions:

######  Definition 2.2.

The matrix $FV^{-1}$ is called the the next generation matrix.

######  Definition 2.3.

The parameter $R_{0}=\rho(FV^{-1})$ is called the basic reproduction number.

In Def. 2.3, $\rho(\cdot)$ denotes the spectral radius of a matrix.

##  3 Evaluating the Effect of a Treatment-type Control

###  3.1 A controlled disease propagation model

We wish to model the effect of treatment-type controls, i.e., the controls
that result in moving individuals from the infected compartments to other
(both infected and healthy) compartments. In doing so, we restrict our
attention to the cases where the control action enters the equations linearly.
In particular, this implies that there is no interference between different
treatment strategies, i.e., the rate at which people are administered to
treatment $i$ does not depend on the respective rate associated with treatment
$j$ for all pairs of $i$ and $j$.

With the above assumptions in mind, we write the controlled population balance
model in the following form:

$\dot{x}=\Phi(x)+\Phi^{u}(x)u,$ (5)

where $u\in\mathbb{R}^{m}_{\geq 0}$ is the vector of non-negative controls and
the components of the control matrix
$\Phi^{u}(x)=\left[\Phi_{ik}^{u}(x)\right]$, $i=1,\dots,n$, $k=1,\dots,m$,
have the following structure:

$\Phi^{u}_{\cdot,k}(x)=\begin{bmatrix}\left(\sum\limits_{\begin{subarray}{c}j\leq
l\\\ j\neq i\end{subarray}}a^{u}_{ij,k}(x)-\sum\limits_{j\neq
i}a^{u}_{ji,k}(x)\right)_{i=1,\dots,l}\\\ \left(\sum\limits_{j\leq
l}a^{u}_{ij,k}(x)\right)_{i=l+1,\dots,n}\end{bmatrix}.$

The control matrix describes the controlled flows from the infected
compartments to both infected and healthy compartments. This may correspond to
isolating infected individuals (corresponds to a transition from an
$I$-compartment to another $I$-compartment) or treating an infected individual
(a flow from an $I$-compartment to an $N$-compartment).

The respective controlled flow rates $a^{u}_{ij,k}(x)$ satisfy the conditions
similar to those for the flow rates of the original (uncontrolled) model:

$\displaystyle a^{u}_{ji,k}(x)\geq 0,$ $\displaystyle\forall
x\in\mathbb{R}^{n}_{\geq 0},i,j=1,\dots,n,\,i\neq j,\,k=1,\dots,m.$ (6a)
$\displaystyle\\{x_{i}=0\\}\implies\\{a^{u}_{ji,k}(x)=0\\}$
$\displaystyle\forall i,j=1,\dots,n,\,k=1,\dots,m.$ (6b)

#### Controllability

Note that at a DFE $x^{*}$, the control matrix vanishes identically as there
are no outflows from I-compartments: $\Phi^{u}(x^{*})=0_{[n\times m]}$. This
implies that at a DFE $x^{*}$, the linearized model is effectively
uncontrollable. Thus, one cannot use the linearized model to design a feedback
control law stabilizing the DFE.

###  3.2 Analysis of the controlled system

Due to the difficulties outlined above, we consider a somewhat simpler, but
even more practically relevant situation. Namely, we wish to study the
system’s dynamics when a constant control $u^{*}\in\mathbb{R}^{m}_{\geq 0}$ is
applied. We first formulate the result that extends Lemma 2.1 to include the
controls $u$.

######  Lemma 3.3.

Let $u(t)=u^{*}\in\mathbb{R}^{m}_{\geq 0}$ for all $t\geq 0$. The linearized
model of (5) at a DFE $x^{*}$ has the form

$\dot{x}=\left(A^{I}+A^{C}+\sum_{k=1}^{m}B_{k}u^{*}_{k}\right)x,$ (7)

where $A^{I}$ and $A^{C}$ are defined in Lemma 2.1 and $B_{k}$, $k=1,\dots,m$
are

$B_{k}=\frac{D\Phi_{\cdot,k}^{u}(x)}{Dx}\bigg{|}_{x=x^{*}}=\begin{bmatrix}B_{11,k}&0\\\\[5.0pt]
B_{21,k}&0\end{bmatrix}.$

The matrices $B_{11,k}$ have non-positive diagonal and non-negative off-
diagonal elements. Furthermore, $B_{11,k}$ are weakly column diagonally
dominant.

###### Proof.

See Appendix A. ∎

For the sake of notational simplicity and to comply with the previously
accepted notation, we will write $W_{i}=-B_{11,i}$. Matrices $W_{i}$ belong to
the class of Z-matrices (see Appendix B for details). However, these are not
M-matrices as matrices $W_{i}$ are typically rank-deficient and hence non-
invertible.

We are interested in the stability of the linearized system (7). As in the
uncontrolled case, the structure matrix of (7) is a block-lower triangular
matrix, whose eigenvalues coincide with the eigenvalues of the diagonal
blocks. The lower right block is Hurwitz according to A3. Thus, the stability
of (7) is determined by the eigenvalues of the matrix

$J_{u}=F-V-\sum_{k=1}^{m}W_{k}u_{k}.$

The following lemma gives an important result that will be used in the sequel.

######  Lemma 3.4.

Let $V$ be a strictly column diagonally dominant Z-matrix with positive
diagonal elements. Then for any $u\in\mathbb{R}^{m}_{\geq 0}$ the matrix
$V+\sum_{k=1}^{m}W_{k}u_{k}$ is an M-matrix.

###### Proof.

According to Lemma 3.3, matrices $W_{i}=-B_{11,i}$ have positive diagonal
elements and are weakly (non-strictly) column diagonally dominant. Thus for
any positive $u$ the matrix $V+\sum_{k=1}^{m}W_{k}u_{k}$ is a Z-matrix and is
strictly column diagonally dominant. Then it follows from Theorem 5.9 that
$V+\sum_{k=1}^{m}W_{k}u_{k}$ is a non-singular M-matrix. ∎

Now we are ready to formulate a generalized version of Theorem 2.2

######  Theorem 3.5.

The matrix $F-V-\sum_{k=1}^{m}W_{k}u_{k}$ is Hurwitz if and only if

$\rho\left(F\left(V+\sum_{k=1}^{m}W_{k}u_{k}\right)^{-1}\right)<1.$

We thus define the controlled reproduction number $R_{0}(u)$ as the spectral
radius of the perturbed matrix
$Q(u)=F\left(V+\sum_{k=1}^{m}W_{k}u_{k}\right)^{-1}$:

$R_{0}(u)=\rho\left(F\left(V+\sum_{k=1}^{m}W_{k}u_{k}\right)^{-1}\right).$ (8)

Obviously, we have $R_{0}(0)=R_{0}$.

Let the uncontrollable system be such that $R_{0}>1$. Theorem 3.5 allows for
determining if a given constant control $u^{*}$ suffices to shift the value of
the basic reproduction number in order to make it less than 1. However, in
many cases it is difficult to compute the perturbed reproduction number
$R_{0}(u)$. Thus, we would like to have a result that would tell us if a given
structure of the treatment allows for achieving the stated goal.

The following result provides the required estimation. We first consider the
case $m=1$, i.e., we assume that there is a scalar control $u>0$.

######  Theorem 3.6.

Let $R_{0}$, $x_{0}$ and $y_{0}$ be the spectral radius of $FV^{-1}$ as well
as the right and the left eigenvectors of $FV^{-1}$ corresponding to $R_{0}$,
respectively. Let, furthermore, there be only one eigenvalue of $FV^{-1}$
coinciding with $R_{0}$ and other eigenvalues be strictly less than $R_{0}$ in
absolute value. For sufficiently small $u$, the sign of variation
$R_{0}(u)-R_{0}$ is determined by the sign of

$R_{1}=-y_{0}^{\top}V^{-1}Wx_{0}R_{0}/(y_{0}^{\top}x_{0})$ (9)

###### Proof.

The eigenvalues of $Q(u)=F\left(V+Wu\right)^{-1}$ change continuously with
$u$. Hence, we can write the spectral radius $R_{0}(u)=R_{0}+uR_{1}+O(u^{2})$
and the respective eigenvector as $x(u)=x_{0}+ux_{1}+O(u^{2})$. We thus have

$F(V+uW)^{-1}(x_{0}+ux_{1}+O(u^{2}))=(R_{0}+uR_{1}+O(u^{2}))(x_{0}+ux_{1}+O(u^{2})).$
(10)

First, we note that $\frac{d}{du}(V+uW)^{-1}=-(V+uW)^{-1}W(V+uW)^{-1}$ and at
$u=0$ we have $\frac{d}{du}(V+uW)^{-1}\bigg{|}_{u=0}=-V^{-1}WV^{-1}$.
Differentiating the left and the right sides of (10) and setting $u=0$ we get:

$(FV^{-1}-R_{0}I)x_{1}=(FV^{-1}WV^{-1}+R_{1}I)x_{0}.$ (11)

The matrix $(FV^{-1}-R_{0}I)$ has a zero eigenvalue and the respective left
eigenvector is $y_{0}$, i.e., $y_{0}^{\top}(FV^{-1}-R_{0}I)=0$. Multiplying
both sides of (11) with $y_{0}^{\top}$ and expressing $R_{1}$ we get:

$R_{1}=-y_{0}^{\top}WV^{-1}x_{0}R_{0}/(y_{0}^{\top}x_{0}).$

Finally, we note that for sufficiently small (and positive) $u$ the sign of
the difference $R_{0}(u)-R_{0}$ is determined by the sign of $R_{1}$, whence
the result follows. ∎

That is to say, we can determine if a given treatment is efficient (at least
for small values of $u$) by checking the sign of $R_{1}$. If this sign is
negative then by increasing $u$ we decrease $R_{0}(u)$ and eventually ensure
that it becomes less than 1. Otherwise (if $R_{1}>0$) we conclude that the
treatment program is inadequately formulated and does not lead to a decrease
in the reproduction number. See Sec. 4 for examples.

Following the same procedure, one can write $R_{0}(u)$ for the case when
$u\in\mathbb{R}^{m}$. We have

$R_{0}(u^{*})=R_{0}+\sum_{k=1}^{m}u^{*}_{k}R^{k}_{1}+O(\|u^{*}\|^{2}).$

If the components of $u^{*}$ are sufficiently small, the contribution of each
individual component is determined by the respective term $R^{k}_{1}$, which
is defined as

$R^{k}_{1}=-y_{0}^{\top}W_{k}V^{-1}x_{0}R_{0}/(y_{0}^{\top}x_{0}).$ (12)

We immediately arrive at the following result.

######  Theorem 3.7.

For sufficiently small values of the controls $u_{i}$, the effect of each
control is independent from the values of the remaining controls. The total
change of $R_{0}(u)$ is equal to the sum of individual contributions up to the
high order term: $R_{0}(u)-R_{0}=\sum_{k=1}^{m}u_{k}R^{k}_{1}+O(\|u\|^{2})$.

This result indicates the second potential use of $R_{1}$. Let there be $m$
treatment strategies that aim at decreasing the value of $R_{0}(u)$. Then the
most efficient strategy is the one with the smallest (and necessarily
negative) value of the parameter $R^{k}_{1}$.

#### Discussion

Note that the structure of the expression for $R_{1}$ does resemble that one
for $R_{0}$ and it can be interpreted as follows. First, following [28], we
interpret the $(i,j)$ entry of $V^{-1}$ as the average duration of time an
individual introduced into the $j$-th compartment spends in the $i$-th one,
assuming there is no reinfection and no control. The $(i,i)$ entry of
$(-W_{k})$ is the rate at which infected individuals are removed from the
$i$-th compartment, while the $(i,j)$-th component, $i\neq j$, corresponds to
the rate at which the infected from $i$-th compartment are moved to the $j$-th
compartment. The respective rates are multiplied by the control $u_{k}$.

Hence, the $(i,j)$-th element of $-WV^{-1}$ is the relative treatment-induced
rate of flow (outflow if negative or inflow if positive) from the $i$-th
compartment as applied to the individuals initially introduced into the $j$-th
compartment. Note that $\sum_{i}(-WV^{-1})_{i,j}=0$ if the treatment consists
in redistributing the infected individuals between infected compartments and
is negative if at least a part of infected individuals are removed to the
susceptible compartments as a result of treatment.

To interpret the right eigenvector of $FV^{-1}$, corresponding to $R_{0}$,
recall that $FV^{-1}x_{0}$ is the expected number of new infections produced
by the initial distribution of infected individuals given by $x_{0}$. Hence,
the right eigenvector $x_{0}$ can be interpreted as the worst case
distribution of the infected. In contrast to that, the left eigenvector
$y_{0}$ can be seen as a worst case transmissibility rates, i.e., the
transmissibility rates that result in the maximal infection spread in the
population taking into account the existing structure of the transmission
routes. Note that the eigenvectors are defined up to a positive factor and
hence should be seen as proportions rather than absolute values.

Finally, the expression $-y_{0}^{\top}WV^{-1}x_{0}(y_{0}^{\top}x_{0})^{-1}$
can be interpreted as the total flow of the infection due to the treatment
computed for the worst case scenario. The term $(y_{0}^{\top}x_{0})^{-1}$
serves as a normalizing factor. If the total flow is negative, the treatment
leads to a decrease in infection and to an increase otherwise. Note that
$R_{1}$ is negative if $\sum_{i}(-WV^{-1})_{i,j}<0$ for all $j$.

######  Remark 3.1.

There are two main advantages of considering $R_{1}$, resp., $R_{1}^{k}$ in
contrast to working directly with $R_{0}(u)$ as defined in (8):

  1. 1.

The function $R_{0}(u)$ depends on $u$ in a highly nonlinear way, making it
difficult if ever possible to analyze the impact of the control on the
controlled reproduction number $R_{0}(u)$. In contrast, $R_{1}$ is a well-
defined quantity that does not depend on $u$ thus making it more amenable for
analysis.

  2. 2.

In general, the matrix $\left(V+\sum_{k=1}^{m}W_{k}u_{k}\right)^{-1}F$ has a
complex structure which substantially complicates the problem of finding its
spectral radius in an analytic form. In contrast to that, $R^{k}_{1}$ can be
computed with much less effort. Also, we note that the vectors $x_{0}$ and
$y_{0}$ need to be computed only once. The modeler can then vary the structure
of $\Phi^{u}$ (and, respectively, $W_{k}$) to achieve the required effect in
$R^{k}_{1}$.

######  Remark 3.2.

Note that the preceding results are formulated for sufficiently small values
of $u$. However, this assumption is not as restrictive as it may seem. The
point is that each component of the control $u$ denotes a fraction of the set
of potentially eligible individuals that is enrolled into the respective
treatment during a unit of time. In most practically relevant cases, this
fraction is rather small, ranging from thousandths to hundredths, thus
justifying the assumption.

###  3.3 Further extensions

The described method is substantially based upon the eigenvalue perturbation
theory. Namely, we analyze the behavior of the largest real eigenvalue
(spectral radius) of a next generation matrix under certain structured
perturbations. All the presented results are of local nature, i.e., these are
valid for sufficiently small values of $u$. However, in many cases it is
possible to extend this analysis to larger values of $u$. Below, we consider
several possible scenarios and consider possible options and potential
limitations.

The simplest case is when the rank of the matrix $F$ is equal to 1. This
corresponds to the situation when all immediately infected individuals enter
the same compartment, but the infection itself can be caused by a contact with
different types of infected. In this case, there is only one non-zero positive
eigenvalue coinciding with the spectral radius. The only possible limitation
corresponding to this situation is that the controlled reproduction number can
change its behavior when the control $u$ grows sufficiently large. One could
thus wish to study the behavior of $R_{0}(u)$ in some more detail. To do so
one can compute the second term in the series expansion which we denote by
$R_{2}$. The corresponding expression is presented in Appendix A, Eq. (25).
However, determination of $R_{2}$ is a rather cumbersome procedure involving
computing Moore-Penrose pseudo inverse. On the other hand, it seems that most
cases the sign of $R_{1}$ describes the behavior of $R_{0}(u)$ for arbitrary
large values of $u$ and hence can be used for making global predictions.

When the rank of $F$ is larger than 1, the matrix $(V+uW)^{-1}F$ has in
general more than one nonzero eigenvalue. A model of co-infection and two
types of treatment, one for each type of infection, is a typical example of
such a situation. In this case, there are two nonzero eigenvalues
$\lambda_{1}(u_{1})$ and $\lambda_{2}(u_{2})$ controlled by the respective
treatments $u_{1}$ and $u_{2}$. Suppose that for $u_{1}=u_{2}=0$,
$R_{0}(0)=\lambda_{1}(0)>\lambda_{2}(0)$. Hence, the optimal strategy (in
terms of minimizing $R_{0}(u)$) would be to invest into $u_{1}$. However, as
$\lambda_{1}(u_{1})$ decreases, it will at some point equate with
$\lambda_{2}(u_{2})$. Further investment into $u_{1}$ will lead to the change
of the roles: $\lambda_{2}(u_{2})$ will become greater than
$\lambda_{1}(u_{1})$ and hence will account for $R_{0}(u)$.

In general, the behavior of the eigenvalues can be predicted by analyzing the
structure of $Q(u)=(V+uW)^{-1}F$. We assume here that this structure does not
change with $u$, that is, $Q_{ij}(u)=0,u>0\Rightarrow Q_{ij}(u)=0,\forall u>0$
and $Q_{ij}(u)\neq 0,u>0\Rightarrow Q_{ij}(u)\neq 0,\forall u>0$. In this case
there are two possible situations: the matrix $Q(u)$ can be either irreducible
or reducible. The former implies that there is a simple real eigenvalue
corresponding to the spectral radius and hence to $R_{0}$ (see, e.g., [24,
Sec. 8.3]). If the matrix $Q(u)$ is reducible, one can find a permutation
matrix $P$ such that $\bar{Q}(u)=PQ(u)P^{-1}$ has a block-diagonal form and
the blocks on the main diagonal are irreducible. Since the eigenvalues of the
permuted matrix are determined by the diagonal blocks $\bar{Q}_{i}(u)$ it
suffices to track the spectral radii of $\bar{Q}_{i}(u)$; denote them by
$\bar{\lambda}_{i}(u)$. We have that
$R_{0}(u)=\max_{i}\left(\bar{\lambda}_{i}(u)\right)$. If for some $u^{*}$ it
happens that $\bar{\lambda}_{i}(u^{*})=\bar{\lambda}_{j}(u^{*})$ we have the
situation described in the preceding paragraph.

##  4 Examples

In this section we consider three sufficiently simple but yet non-trivial
epidemiological models that are aimed at illustrating the usefulness of the
proposed approach.

###  4.1 An SI model with acute and chronic stages and a single treatment

#### Model

Consider an SI-model of a disease with two stages (acute and chronic) and a
treatment $u$. This general model structure with sequential infection
compartments containing infected persons with different levels of
contagiousness describes the evolution of such diseases as HIV [25] and
syphilis [13] that have variable levels of contagiousness over the course of
infection. For the sake of illustration we explicitly show the $\Phi^{I}$,
$\Phi^{C}$, and $\Phi^{u}$ components in (13).

$\frac{d}{dt}\begin{bmatrix}I_{A}\\\ I_{C}\\\ T\\\
S\end{bmatrix}=\underbrace{\begin{bmatrix}\alpha(X)S\\\ 0\\\ 0\\\
0\end{bmatrix}}_{\Phi^{I}(x)}+\underbrace{\begin{bmatrix}-\beta
I_{A}-\mu_{A}I_{A}\\\ \beta I_{A}+\delta T-\mu_{C}I_{C}\\\ -\delta
T-\mu_{T}T\\\
w-\alpha(X)S-\mu_{S}S\end{bmatrix}}_{\Phi^{C}(x)}+\underbrace{\begin{bmatrix}0\\\
-I_{C}\\\ I_{C}\\\ 0\end{bmatrix}}_{\Phi^{u}(x)}u$ (13)

where $\alpha(X)=\dfrac{\alpha_{A}I_{A}+\alpha_{C}I_{C}+\alpha_{T}T}{N}$,
$X=[I_{A},I_{C},T,S]$, and $N=I_{A}+I_{C}+T+S$. The states correspond to the
number of acutely infected ($I_{A}$), chronically infected ($I_{C}$), treated
($T$), and susceptible ($S$) individuals. We assume that the first three
compartments are infectious with different transmission probabilities:
$\alpha_{A}$, $\alpha_{C}$, and $\alpha_{T}$. Furthermore, $\beta$ is the
inverse duration of the acute phase, $\mu_{S}$, $\mu_{A}$, $\mu_{C}$ and
$\mu_{T}$ are the mortality rates for susceptible, non-treated in acute and
chronic phases and treated infected, with
$\mu_{C}>\\{\mu_{A},\mu_{T}\\}>\mu_{S}$, and $\delta$ is the rate at which the
treatment fails.

#### Analysis

The disease-free equilibrium is unique and given by $x^{*}=[0,0,0,w/\mu_{S}]$.
The (uncontrolled) value of $R_{0}$ is computed as the spectral radius of
$FV^{-1}$ and is equal to

$R_{0}=\frac{\alpha_{C}\beta+\alpha_{A}\mu_{C}}{\left(\beta+\mu_{A}\right)\mu_{C}}$

Note that the DFE $x^{*}$ is unstable if
$\alpha_{C}\beta+\alpha_{A}\mu_{C}>(\beta+\mu_{S})\mu_{C}$. We wish to
evaluate if the proposed treatment strategy is efficient. To do so we compute
$R_{1}$ according to (9):

$R_{1}=-\frac{\beta\left(\alpha_{C}\mu_{T}-\alpha_{T}\mu_{C}\right)}{\mu_{C}^{2}\left(\beta+\mu_{A}\right)\left(\delta+\mu_{T}\right)}$
(14)

This expression is negative if $\alpha_{C}\mu>\alpha_{T}\mu_{C}$. This can be
written as

$\alpha_{C}\mu_{C}^{-1}>\alpha_{T}\mu_{T}^{-1},$ (15)

that is the treatment turns out to be efficient if the product (transmission
probability times average residence time in the respective compartment) is
lower for the treated individuals as compared with chronically infected. In
the following, we will refer to this product as the cumulative
transmissibility of the infection induced by a given compartment. The
condition (15) can be interpreted as follows: If treated people live longer
than untreated (but infected), but their infectivity does not reduce enough to
compensate for this, such treatment will contribute to the propagation of the
disease.

For the considered model it is possible to compute the controlled reproduction
number $R_{0}(u)$:

$R_{0}(u)=\frac{\left(\alpha_{C}\beta+\alpha_{A}\mu_{C}\right)\left(\delta+\mu_{T}\right)+\left(\alpha_{T}\beta+\alpha_{A}\mu_{T}\right)u}{(\beta+\mu_{A})(\delta\mu_{C}+\mu_{C}\mu_{T}+\mu_{T}u)}.$

One may check that $R_{1}=\dfrac{\partial R_{0}(u)}{\partial
u}\bigg{|}_{u=0}$. However, we are interested in the behavior of $R_{0}(u)$
for (relatively) large values of $u$. To do so, we compute the difference

$R_{0}(u)-R_{0}=-\frac{\beta
u\left(\alpha_{C}\mu_{T}-\alpha_{T}\mu_{C}\right)}{\mu_{C}\left(\beta+\mu_{A}\right)\left(\delta\mu_{C}+\mu_{C}\mu_{T}+\mu_{T}u\right)}.$
(16)

One can readily observe that the difference is negative for any $u>0$ when
(15) holds.

Figure 1 shows the relative error in the $R_{1}$ approximation to $R_{0}(u)$
for a set of random parameters as a function of $R_{0}$. We obtained random
parameter sets for a given value of $R_{0}$ optimizing from a random starting
point for $\alpha_{A}$, $\alpha_{C}$, and $\beta$ assuming
$\mu_{A}=\mu_{C}=120^{-1}$, $\mu_{T}=360^{-1}$, and $\alpha_{T}=0$. We
rejected parameter sets such that there was no feasible control capable of
reducing $R_{0}(u)$ to 1. The relative error is defined as the difference in
the exact value of $u^{*}$ such that $R_{0}(u^{*})=1$ and
$u^{**}=\frac{R_{0}-1}{R_{1}}$ normalized by $u^{*}$. The approximate value of
the control is always less than the true value and the approximation becomes
worse as the value of $R_{0}$ increases, although the accuracy of the
approximation is dependent on the specific parameter values. However, if
$R_{0}$ is small the approximation gives a reasonable indicator of the level
of control required to bring the epidemic to the threshold level.

Figure 1: Relative error in approximate reduction in $R_{0}$ due to treatment.
This figure shows the multiplicative-scale error defined as
$\frac{u^{*}-u^{**}}{u^{*}}$ where $u^{*}$ such that $R_{0}(u^{*})=1$ and
$u^{**}=\frac{1-R_{0}}{R_{1}}$ for a set of randomly selected parameter
values. Parameters with a given $R_{0}$ between 1.1 and 4 were selected at
random as described in the text. The blue line is a local polynomial
regression showing the general trend.

###  4.2 An SEIR model with an asymptomatic stage and treatment

#### Model

Consider the model shown in Fig. 2. This model describes the transmission
dynamics of an asymptomatic (sub-clinical) infection. Individuals with
asymptomatic infection do not develop the respective symptoms, but are
infectious and contribute to the distribution of the disease. Asymptomatic
infection has been shown to exist for many diseases, including, e.g., herpes
[29], gonorrhea [20], measles [1], and common cold.

Figure 2: SEIR model.

Here, we have the following compartments: $E$ – exposed, $I$ – infected, $A$ –
asymptomatically infected, $T$ – treated. $\Lambda$ is the inflow into the
population, $\eta$ is the transmission rate, $p_{1}$ is the probability that
the exposition will lead to the typical clinical course and $p_{2}$ is the
probability that the failed treatment will result in the typical clinical
course. Furthermore, $x$ is the rate at which the treatment is failed or
canceled, and $\gamma_{I}$, $\gamma_{A}$ and $\gamma_{T}$ are the inverses of
the mean residence time in respective compartments. We assume that only the
symptomatic infected, i.e., $I$, are treated. However, the disease can be
transmitted both by symptomatic and asymptomatic infected though with
different transmission rates $\alpha_{I}$ and $\alpha_{A}$. People under
treatment are assumed to be non-infectious.

Its behavior is described by the following DEs:

$\begin{cases}\dot{E}=\alpha(X)S-(\eta+\mu)E\\\
\dot{A}=-(\gamma_{A}+\mu)A+\eta(1-p_{1})E+x(1-p_{2})T\\\ \dot{I}=\eta
p_{1}E-\left(\gamma_{I}+\mu\right)I-uI+p_{2}xT\\\
\dot{T}=uI-(\gamma_{T}+\mu+x)T\\\ \dot{S}=\Lambda-\mu S-\alpha(X)S\\\
\dot{R}=\gamma_{A}A+\gamma_{I}I+\gamma_{T}T-\mu R,\end{cases}$

where $\alpha(X)=\dfrac{\alpha_{A}A+\alpha_{I}I}{N}$.

#### Analysis

The basic reproduction number $R_{0}$ is found using the NGM approach to be

$R_{0}=\frac{\eta}{\eta+\mu}\left[\frac{\alpha_{A}}{\gamma_{A}+\mu}(1-p_{1})+\frac{\alpha_{I}}{\gamma_{I}+\mu}p_{1}\right],$

where the expression in square brackets is the expected value of the total
cumulative transmissibility of the infection and the quantity in front of the
brackets is the fraction of the infected that leave $E$ toward either $I$ or
$A$.

We compute $R_{1}$

$R_{1}=-\dfrac{\eta
p_{1}\big{[}\alpha_{I}(\gamma_{A}+\mu)(\gamma_{T}+\mu)+\left(\alpha_{I}(\gamma_{A}+\mu)-\alpha_{A}(\gamma_{I}+\mu)\right)(1-p_{2})x\big{]}}{\left(\eta+\mu\right)\left(\gamma_{A}+\mu\right){\left(\gamma_{I}+\mu\right)}^{2}\left(\gamma_{T}+\mu+x\right)}$
(17)

which turns out to be negative if the numerator is positive (note that the
denominator is always positive).

Upon some algebraic manipulations, this condition can be written as

$\frac{\alpha_{I}}{\gamma_{I}+\mu}>\frac{x(1-p_{2})}{\gamma_{T}+\mu}\left[\frac{\alpha_{A}}{\gamma_{A}+\mu}-\frac{\alpha_{I}}{\gamma_{I}+\mu}\right].$
(18)

Define the cumulative transmissibilities of the $A$ and $I$ stages as
$\zeta_{A}=\frac{\alpha_{A}}{\gamma_{A}+\mu}$ and
$\zeta_{I}=\frac{\alpha_{I}}{\gamma_{I}+\mu}$. With this, the condition (18)
can be rewritten as

$\frac{\zeta_{A}-\zeta_{I}}{\zeta_{A}}\cdot x(1-p_{2})<\gamma_{T}+\mu,$

provided $\zeta_{A}>\zeta_{I}$ (otherwise the condition (18) holds trivially
as the r.h.s. of (18) turns to be negative). The first term on the left side
is the relative difference in cumulative transmissibility and the second one
is the flow from $T$ to $A$. Finally, the expression on the right side
describes the flow from $T$ to $R$ or to outside the system. This result
underscores the importance of the treatment efficiency for the successful
eradication of the disease.

We proceed by considering the total change of $R_{0}$ due to the control $u$

$R_{0}(u)-R_{0}=\\\ -\dfrac{\eta
p_{1}u\left[\alpha_{I}(\gamma_{A}+\mu)(\gamma_{T}+\mu)+\left(\alpha_{I}(\gamma_{A}+\mu)-\alpha_{A}(\gamma_{I}+\mu)\right)(1-p_{2})x\right]}{(\eta+\mu)(\gamma_{A}+\mu)(\gamma_{I}+\mu)\left[(\gamma_{I}+\mu)(\gamma_{T}+\mu+x)+(\gamma_{T}+\mu+x(1-p_{2}))u\right]}.$
(19)

The denominator of (19) is always positive so, the condition for the treatment
to be efficient (that is the condition for the difference to be negative) is
that the numerator is positive. But the numerator of (19) coincides with the
numerator of $R_{1}$, (17), up to the control $u$ which is positive. Hence we
conclude that satisfaction of the condition (18) guarantees decrease in
$R_{0}$ for any positive value of $u$, i.e., gives a global condition.

###  4.3 An SI model with high- and low-risk groups and two treatments

#### Model

Consider a simplified model of HIV transmission dynamics with two controls
corresponding to treating infected individuals from the high- and the low-risk
groups, denoted by $u_{T}$ and $u_{L}$. For the details on the derivation of
the model see [18]. The transmission dynamics is described by the following
set of ODEs

$\begin{cases}\dot{I}_{H}=-I_{H}\mu-
I_{H}f_{L}\omega+I_{L}f_{H}\omega+\phi_{H}(X)S_{H}-I_{H}u_{H}-\mu_{I}I_{H}\\\
\dot{I}_{L}=-I_{L}\mu+I_{H}f_{L}\omega-
I_{L}f_{H}\omega+\phi_{L}(X)S_{L}-I_{L}u_{L}-\mu_{I}I_{L}\\\
\dot{T}_{H}=T_{L}f_{H}\omega-T_{H}\mu-T_{H}f_{L}\omega+I_{H}u_{H}\\\
\dot{T}_{L}=T_{H}f_{L}\omega-T_{L}\mu-T_{L}f_{H}\omega+I_{L}u_{L}\\\
\dot{S}_{H}=\alpha_{H}-S_{H}\mu-
S_{H}f_{L}\omega+S_{L}f_{H}\omega-\phi_{H}(X)S_{H}\\\
\dot{S}_{L}=\alpha_{L}-S_{L}\mu+S_{H}f_{L}\omega-
S_{L}f_{H}\omega-\phi_{L}(X)S_{L},\end{cases}$ (20)

where $S_{\\{H,L\\}}$, $I_{\\{H,L\\}}$, and $T_{\\{H,L\\}}$ are the
susceptible, infected, and treated. The subscript denotes the behavioral
pattern of the respective group: there are a (H)igh and a (L)ow risk group.
Further,
$\phi_{H}(X)=\beta\lambda_{H}\dfrac{\lambda_{H}I_{H}+\lambda_{L}I_{L}}{\lambda_{H}N_{H}+\lambda_{L}N_{L}}$
and
$\phi_{L}(X)=\beta\lambda_{L}\dfrac{\lambda_{H}I_{H}+\lambda_{L}I_{L}}{\lambda_{H}N_{H}+\lambda_{L}N_{L}}$
are the per-capita transmission rates. $\mu$ is the mortality rate and
$\mu_{I}$ is the disease-induced mortality; $\rho_{H}=f_{H}\omega$ and
$\rho_{L}=(1-f_{H})\omega=f_{L}\omega$ are the transition rates between high-
and low-risk groups with $\omega$ denoting the volatility coefficient;
$\lambda_{H}$ and $\lambda_{L}$ are the contact rates; $\beta$ is the
infection transmissibility. Finally, the inflow rates are $\alpha_{H}=f_{H}\mu
N$ and $\alpha_{L}=f_{L}\mu N$.

For the considered model, the disease free equilibrium is
$[I_{H}^{*},\,I_{L}^{*},\,T_{H}^{*},\,T_{L}^{*},\,S_{H}^{*},\,S_{L}^{*}]=[0,\,0,\,0,\,0,\,f_{H}N,\,f_{L}N]$.
Note that $f_{H}$ and $f_{L}$ are the fractions of the respective (high- or
low-risk) population at the DFE. Using the next generation matrix method we
can compute $R_{0}$:

#### Analysis

The basic reproduction number is computed using the NGM method:

$R_{0}=\frac{\beta}{\left(\mu_{I}+\mu+\omega\right)}C_{V}(\lambda)+\frac{\beta}{\left(\mu_{I}+\mu\right)}[\lambda]$

where $[\lambda]=f_{L}\lambda_{L}+f_{H}\lambda_{H}$ and
$[\lambda^{2}]=\lambda_{L}^{2}f_{L}+\lambda_{H}^{2}f_{H}$ are the first and
the second moments of the contact rate at the DFE,
$\operatorname{Var}(\lambda)$ is the variance of the contact rate, and $C_{V}$
is the coefficient of variance defined as
$C_{V}(\lambda)=\operatorname{Var}(\lambda)\bigm{/}[\lambda]$. Since there are
two controls, we would like to compare their contributions in order to decide
which one should be invested into. To compute the corresponding components
$R_{1}^{H}$ and $R_{1}^{L}$ we use (12) to get

$R_{1}^{H}=-\frac{\beta
f_{H}{\left(\lambda_{H}(\mu_{I}+\mu)+[\lambda]\omega\right)}^{2}}{\left(\mu_{I}+\mu\right)^{2}[\lambda]{\left(\mu_{I}+\mu+\omega\right)}^{2}}$

$R_{1}^{L}=-\frac{\beta
f_{L}\left(\lambda_{L}(\mu_{I}+\mu)+[\lambda]\omega\right)^{2}}{{\left(\mu_{I}+\mu\right)}^{2}[\lambda]{\left(\mu_{I}+\mu+\omega\right)}^{2}}$

The first observation is that both $R_{1}^{H}$ and $R_{1}^{L}$ are negative
thus, they contribute to reducing $R_{0}(u)$ for any choice of parameters.
After some algebraic manipulations we find that $u_{H}$ is more efficient than
$u_{L}$ if

$f_{H}{\left(\frac{\lambda_{H}}{[\lambda]}+\frac{\omega}{\mu}\right)}^{2}>f_{L}\left(\frac{\lambda_{L}}{[\lambda]}+\frac{\omega}{\mu}\right)^{2}$

Figure 3 shows the value of $\lambda_{H}$ such that the two terms in the
inequality above are equal. To account for the variability in the duration of
the high and low-risk periods as a function of $f_{H}$, we introduce the
normalized volatility coefficient

$\omega^{*}=\frac{\omega\mu^{-1}}{{f_{H}}^{-1}+(1-f_{H})^{-1}},$

which is the number of full high and low-risk episodes that can be contained
in a typical infectious period. This plot shows that for fixed value of
$f_{H}$, behavioral volatility makes high-risk intervention more plausible
(i.e. the high-risk population does not have to be extremely high-risk to make
a targeted intervention efficient). This simple analysis gives us a clear
theoretical prescription for when to focus on high-risk group based on
measurable aspects of the transmission system.

Figure 3: Minimum value of the high-risk contact rate, $\lambda_{H}$, for
which the high-risk intervention is preferred. The color shows the value of
$\lambda_{H}$ for which the terms in the above inequality are equal. For any
value of $\lambda_{H}$ higher than the plotted value the high-risk
intervention is preferred. The gray color indicates that the value is either
above $10^{4}$ or there is no value such that the terms are equal. The x-axis
is the normalized volatility coefficient $\omega^{*}$. The y-axis shows the
fraction of the population that is high-risk in the absence of disease. The
remaining parameters are $\lambda_{L}=1$ and $\mu_{I}=0$.

##  5 Discussion

The results presented in this paper can be applied to a wide class of
epidemiological models as long as their dynamics can be described by a
compartmental system of form (5). Our approach builds upon and further
develops the next generation matrix method in that it allows one to estimate
the influence of the treatment-type control(s) on the basic reproduction
number $R_{0}$ which defines the ultimate condition for eventual elimination
of a disease. The more a given control reduces $R_{0}$, the closer the system
is to elimination and the more effective future interventions will be.
Furthermore, it may turn out that in complex models an intervention could
unintentionally make things worse for certain populations. The basic premise
of medicine to do no harm applies to public health as well. However, the
complex, non-linear dynamics of transmission limit the ability of our
intuitions to predict the effects of an intervention. Likewise, measurement of
the effects of interventions are often very noisy and can have long time lags.
Both weak measurability of outcomes and hard to predict dynamics highlight the
need for stronger theoretical guarantees that an intervention will not cause
population-level harm. A possible extension to this work could include
consideration of complex models of how risk behavior changes in response to
changing prevalence and incidence of disease.

All our results are of local nature, i.e., these are valid for sufficiently
small values of $u$. However, in many cases it is possible to extend this
analysis to larger values of $u$ as was illustrated by examples It should also
be noted that the proposed approach does inherit all the limitations
associated with the NGM method. For instance, it provides only a local
stability condition and does not allow to make a conclusion about the behavior
of the system under large deviations.

## Acknowledgment

Research presented in this article was supported by the Laboratory Directed
Research and Development program of Los Alamos National Laboratory under
project number 20180612ECR.

## Appendix A. Proofs and computations

#### Proof of Lemma 3.3

One can readily observe that $\dfrac{\partial\Phi^{u}_{i,k}(x)}{\partial
x_{j}}=0$ for all $j>l$ as $\Phi^{u}_{i,k}(x^{*})$ turns to zero identically
for any $x^{*}=\begin{bmatrix}\mathbf{0}_{l}&x^{C}\end{bmatrix}^{\top}$. We
thus consider the partial derivatives of $\Phi^{u}_{i,k}(x^{*})$ w.r.t.
$x_{j}$ for $j=1,\dots,l$.

Let $i,j\in\\{1,\dots,l\\}$. We have

$\displaystyle\dfrac{\partial\Phi^{u}_{i,k}(x)}{\partial x_{j}}=$
$\displaystyle{}\lim_{\delta x_{j}\rightarrow
0}\dfrac{\sum\limits_{\begin{subarray}{c}q\leq l\\\ q\neq
i\end{subarray}}a^{u}_{iq,k}(0,\dots,\delta x_{j},\dots,0)-\sum\limits_{q\neq
i}a^{u}_{qi,k}(0,\dots,\delta x_{j},\dots,0)}{\delta x_{j}}$ $\displaystyle=$
$\displaystyle{}\begin{cases}\lim\limits_{\delta x_{j}\rightarrow
0}\dfrac{-\sum\limits_{q\neq i}a^{u}_{qi,k}(0,\dots,\delta
x_{j},\dots,0)}{\delta x_{j}}\leq 0,&i=j\\\\[10.0pt] \lim\limits_{\delta
x_{j}\rightarrow 0}\dfrac{a^{u}_{ij,k}(0,\dots,\delta x_{j},\dots,0)}{\delta
x_{j}}\geq 0&i\neq j\end{cases}$

For $i\in\\{l+1,\dots,n\\}$ and $j\in\\{1,\dots,l\\}$, the partial derivatives
are

$\displaystyle\dfrac{\partial\Phi^{u}_{i,k}(x)}{\partial x_{j}}=$
$\displaystyle{}\lim_{\delta x_{j}\rightarrow
0}\dfrac{a^{u}_{ij,k}(0,\dots,\delta x_{j},\dots,0)}{\delta x_{j}}\geq 0,$

which yields the required sign structure. Summation over $j$ gives
$\sum_{j=1}^{n}\dfrac{\partial\Phi^{u}_{i,k}(x)}{\partial x_{j}}=0$ for all
$i=1,\dots,l$ and hence, we have (note that the summation is performed only
for $j\leq l$)

$\left|\dfrac{\partial\Phi^{u}_{i,k}(x)}{\partial
x_{j}}\right|\geq\sum\limits_{\begin{subarray}{c}j\leq l\\\ j\neq
i\end{subarray}}\left|\dfrac{\partial\Phi^{u}_{i,k}(x)}{\partial
x_{j}}\right|.$ (21)

This implies that the matrix
$B_{11,k}=\left[\dfrac{\partial\Phi^{u}_{i,k}(x)}{\partial
x_{j}}\right]_{\begin{subarray}{c}i=1,\dots,l\\\ j=1,\dots,l\end{subarray}}$
is weakly column diagonally dominant as the inequality in (21) is not strict.
$\Box$

#### Computation of the second order term in the expansion of $R_{0}(u)$

Let $A(u)$ be a matrix depending on $u$, $r_{0}$ be the simple eigenvalue
equal to the spectral radius of $A(0)$, and $w_{0}$ and $v_{0}$ be the left
and the right eigenvectors corresponding to $r_{0}$. We expand the perturbed
eigenvalue $r(u)$ and the corresponding right eigenvector $x(u)$ in a Taylor
series and keep the terms up to the second order:
$r(u)=r_{0}+r_{1}u+r_{2}u^{2}+O(u^{3})$ and
$x(u)=v_{0}+v_{1}u+v_{2}u^{2}+O(u^{3})$. Thus we have

$A(u)(v_{0}+v_{1}u+v_{2}u^{2}+O(u^{3}))=(v_{0}+v_{1}u+v_{2}u^{2}+O(u^{3}))(r_{0}+r_{1}u+r_{2}u^{2}+O(u^{3}))$
(22)

Differentiating (22) w.r.t. $u$ and evaluating at $u=0$ we get

$(Ir_{0}-A(0))v_{1}=(A^{\prime}(0)-Ir_{1})v_{0},$ (23)

whence the expression for $r_{1}$ can be obtained:
$r_{1}=w_{0}^{\top}A^{\prime}(0)v_{0}(w_{0}^{\top}v_{0})^{-1}$ (cf. the proof
of Thm. 3.6). Substituting $r_{1}$ back to (23) one gets an expression that
can be used to determine $v_{1}$ (see [22, Chap. 8] for details):

$v_{1}=(Ir_{0}-A(0))^{\dagger}\left(I-\frac{v_{0}w_{0}^{\top}}{w_{0}^{\top}v_{0}}\right)A^{\prime}(0)v_{0},$
(24)

where $()^{\dagger}$ is the Moor-Penrose inverse operator.

To compute the second term in the expansion of $r(u)$ we differentiate (22)
twice w.r.t. $u$ and evaluate at $u=0$ to get

$\frac{1}{2}A^{\prime\prime}(0)v_{0}+A^{\prime}(0)v_{1}-v_{1}r_{1}-v_{0}r_{2}=(Ir_{0}-A(0))v_{2}$

Multiplying from the left by $w_{0}^{\top}$ and substituting the previously
obtained expressions for $r_{1}$ and $v_{1}$ we arrive after some computations
to the final expression for $r_{2}$:

$r_{2}=\left(w_{0}^{\top}v_{0}\right)^{-1}w_{0}^{\top}\left[\frac{1}{2}A^{\prime\prime}(0)+A^{\prime}(0)P_{0}(Ir_{0}-A(0))^{\dagger}P_{0}A^{\prime}(0)\right]v_{0},$
(25)

where $P_{0}=I-\frac{v_{0}w_{0}^{\top}}{w_{0}^{\top}v_{0}}$ is the oblique
projection operator. Finally, we recall that $A(u)=(V+uW)^{-1}F$, whence
$A^{\prime}(0)=V^{-1}WV^{-1}F$ and $A^{\prime\prime}(0)=-2V^{-2}WV^{-2}F$.

## Appendix B. Special classes of matrices and their properties.

######  Definition 5.1.

A matrix $B\in\mathbb{R}^{n\times n}$ is non-negative, denoted by $B\succeq
0$, if $b_{ij}\geq 0$ for all $i,j=1,\dots,n$.

######  Definition 5.2.

A non-singular matrix $A$ is said to be inverse positive if it satisfies
$\mathbb{R}^{n}_{\geq 0}\subseteq A\mathbb{R}^{n}_{\geq 0}$, which is
equivalent to $A^{-1}\succeq 0$, i.e., $A^{-1}\mathbb{R}^{n}_{\geq
0}\subseteq\mathbb{R}^{n}_{\geq 0}$.

######  Theorem 5.8.

Let $A=\alpha I-B$, where $\alpha>0$ and $B\succeq 0$. Then the following
statements are equivalent:

  1. 1.

The matrix $A$ is inverse positive,

  2. 2.

The spectral radius of $B$ is strictly smaller than $\alpha$,

  3. 3.

The matrix $A$ is positive stable, i.e., if $\lambda$ is an eigenvalue of $A$,
then $\Re(\lambda)>0$.

######  Definition 5.3.

A matrix $A=\alpha I-B$ satisfying any of the properties of Thm. 5.8 is said
to be an $M$-matrix.

We will occasionally write $\mathbf{M}$ to denote the class of all $n$-by-$n$
non-singular $M$-matrices. The preceding results can be generalized in the
following way. Let us define the class of $Z$-matrices as
$\mathbf{Z}=\\{A\in\mathbb{R}^{n\times n}|a_{ij}\leq 0,i\neq j\\}$. The
following theorem gives a number of conditions which guarantee that a given
Z-matrix is a non-singular M-matrix. For a complete list see [3].

######  Theorem 5.9.

Let $A\in\mathbf{Z}$. Any of the following conditions implies
$A\in\mathbf{M}$.

  1. 1.

$A$ is inverse-positive.

  2. 2.

$A$ is positive stable.

  3. 3.

$A$ has all positive diagonal elements and is strictly row diagonally dominant
(d.d.), i.e.,

$a_{ii}>\sum_{i\neq j}|a_{ij}|,\quad i=1,\dots,n.$

  4. 4.

$A$ has all positive diagonal elements and is strictly column d.d., i.e.,

$a_{ii}>\sum_{i\neq j}|a_{ji}|,\quad i=1,\dots,n.$

###### Proof.

We will prove only the last item as the remaining ones are covered in [3].

Let $A\in\mathbf{Z}$ and $A$ be strictly column d.d., then
$A^{\top}\in\mathbf{Z}$ and is strictly row d.d. This implies
$A^{\top}\in\mathbf{M}$. Since the spectrum of $A$ coincides with that of
$A^{\top}$, the positive stability property holds for $A$ and hence
$A\in\mathbf{M}$. ∎

## References

  * [1] B Anlar, A Ayhan, H Hotta, M Itoh, D Engin, S Barun, and Ö Köseoglu.  Measles virus RNA in tonsils of asymptomatic children.  Journal of paediatrics and child health, 38(4):424–425, 2002.
  * [2] Rebecca F. Baggaley, Neil M. Ferguson, and Geoff P. Garnett.  The epidemiological impact of antiretroviral use predicted by mathematical models: a review.  Emerging Themes in Epidemiology, 2(9), 2005.
  * [3] Abraham Berman and Robert J. Plemmons.  Nonnegative matrices in the mathematical sciences.  SIAM, 1994.
  * [4] S. Blower, E. J. Schwartz, and J. Mills.  Forecasting the future of HIV epidemics: the impact of antiretroviral therapies & imperfect vaccines.  AIDS reviews, 5(2):113–125, 2003.
  * [5] M. C. J. Bootsma and N. M. Ferguson.  The effect of public health measures on the 1918 influenza pandemic in U.S. cities.  Proceedings of the National Academy of Sciences, 104(18):7588–7593, 2007.
  * [6] Francis J Bowden and Geoffrey P Garnett.  Trichomonas vaginalis epidemiology: parameterising and analysing a model of treatment interventions.  Sexually Transmitted Infections, 76(4):248–256, 2000.
  * [7] Carlos Castillo-Chavez and Zhilan Feng.  To treat or not to treat: the case of tuberculosis.  Journal of mathematical biology, 35(6):629–656, 1997.
  * [8] Carlos Castillo-Chavez, Zhilan Feng, and Wenzhang Huang.  On the computation of $R_{0}$ and its role on global stability.  In D. Bies, S. Blower, P. van den Driessche, D. Kirschner, and A.-A. Yakubu, editors, Mathematical approaches for emerging and reemerging infectious diseases: an introduction, volume 125 of The IMA Volumes in Mathematics and its Applications, pages 229–250. Springer, 2002.
  * [9] Earl A. Coddington and Norman Levinson.  Theory of ordinary differential equations.  Tata McGraw-Hill Education, 1955.
  * [10] Odo Diekmann, Hans Heesterbeek, and Tom Britton.  Mathematical tools for understanding infectious disease dynamics.  Princeton University Press, 2012.
  * [11] Odo Diekmann and Johan Andre Peter Heesterbeek.  Mathematical epidemiology of infectious diseases: model building, analysis and interpretation, volume 5 of Mathematical and computational biology.  John Wiley & Sons, 2000.
  * [12] Odo Diekmann, Johan Andre Peter Heesterbeek, and Johan A.J. Metz.  On the definition and the computation of the basic reproduction ratio $r_{0}$ in models for infectious diseases in heterogeneous populations.  Journal of mathematical biology, 28(4):365–382, 1990.
  * [13] Patrick French.  Syphilis.  BMJ : British Medical Journal, 334(7585):143–147, January 2007\.
  * [14] Dmitry Gromov, Ingo Bulla, Oana Silvia Serea, and Ethan O. Romero-Severson.  Numerical optimal control for HIV prevention with dynamic budget allocation.  Mathematical Medicine and Biology: A Journal of the IMA, page dqx015, 2017.
  * [15] J.A.P. Heesterbeek.  A brief history of $R_{0}$ and a recipe for its calculation.  Acta Biotheoretica, 50(3):189–204, Sep 2002.
  * [16] J.A.P. Heesterbeek and K. Dietz.  The concept of $R_{0}$ in epidemic theory.  Statistica Nederlandica, 50(1):89–110, 1996.
  * [17] J.M. Heffernan, R.J. Smith, and L.M. Wahl.  Perspectives on the basic reproductive ratio.  Journal of the Royal Society Interface, 2(4):281–293, 2005.
  * [18] Christopher J Henry and James S Koopman.  Strong influence of behavioral dynamics on the ability of testing and treating HIV to stop transmission.  Scientific reports, 5:9467, 2015.
  * [19] Roger A. Horn and Charles R. Johnson.  Topics in matrix analysis.  Cambridge University Presss, Cambridge, 1991.
  * [20] Eline L Korenromp, Mondastri K Sudaryo, Sake J de Vlas, Ronald H Gray, Nelson K Sewankambo, David Serwadda, Maria J Wawer, and J Dik F Habbema.  What proportion of episodes of gonorrhoea and chlamydia becomes symptomatic?  International Journal of STD & AIDS, 13(2):91–101, 2002.
  * [21] Eline L Korenromp, Carina Van Vliet, Heiner Grosskurth, Awene Gavyole, Catharina PB Van der Ploeg, Lieve Fransen, Richard J Hayes, and J Dik F Habbema.  Model-based evaluation of single-round mass treatment of sexually transmitted diseases for HIV control in a rural African population.  Aids, 14(5):573–593, 2000.
  * [22] Jan R. Magnus and Heinz Neudecker.  Matrix differential calculus with applications in statistics and econometrics.  Wiley, 2nd rev. edition, 1999.
  * [23] Natasha K Martin, Peter Vickerman, Jason Grebely, Margaret Hellard, Sharon J Hutchinson, Viviane D Lima, Graham R Foster, John F Dillon, David J Goldberg, Gregory J Dore, et al.  Hepatitis C virus treatment for prevention among people who inject drugs: modeling treatment scale-up in the age of direct-acting antivirals.  Hepatology, 58(5):1598–1609, 2013.
  * [24] Carl D Meyer.  Matrix analysis and applied linear algebra.  SIAM, 2000.
  * [25] World Health Organization et al.  Interim who clinical staging of hvi/aids and hiv/aids case definitions for surveillance: African region.  Technical report, Geneva: World Health Organization, 2005.
  * [26] Oluwaseun Sharomi, Chandra N. Podder, Abba B. Gumel, and Baojun Song.  Mathematical analysis of the transmission dynamics of HIV/TB coinfection in the presence of treatment.  Mathematical Biosciences & Engineering, 5(1):145–174, 2008.
  * [27] Cristiana J. Silva and Delfim F. M. Torres.  A TB-HIV/AIDS coinfection model and optimal control treatment.  Discrete & Continuous Dynamical Systems - A, 35(9):4639–4663, 2015\.
  * [28] Pauline Van den Driessche and James Watmough.  Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission.  Mathematical biosciences, 180(1):29–48, 2002.
  * [29] Anna Wald, Judith Zeh, Stacy Selke, Rhoda L. Ashley, and Lawrence Corey.  Virologic characteristics of subclinical and symptomatic genital herpes infections.  New England Journal of Medicine, 333(12):770–775, 1995.