# Impact of a cost functional on the optimal control and the cost-
effectiveness: control of a spreading infection as a case study

Fernando Saldaña  José Ariel Camacho-Gutiérrez<EMAIL_ADDRESS>Andrei Korobeinikov  Instituto de Matemáticas, Campus Juriquilla, 76230,
Universidad Nacional Autónoma de México, Quéretaro, Mexico  Facultad de
Ciencias, Universidad Autónoma de Baja California, 22860 Baja California,
Mexico  School of Mathematics and Information Science, Shaanxi Normal
University, Xián, China

###### Abstract

In applications of the optimal control theory to problems in medicine and
biology, the dependency of the objective functional on the control itself is
often a matter of controversy. In this paper, we explore the impact of the
dependency using reasonably simple _SIR_ and _SEIRS_ epidemic models. To
qualitatively compare the outcomes for different objective functionals, we
apply the cost-effectiveness analysis. Our result shows that, at least for the
comparatively inexpensive controls, the variation of the power at the controls
in a biologically feasible range does not significantly affect the forms of
the optimal controls and the corresponding optimal state solutions. Moreover,
the costs and effectiveness are affected even less. At the same time, the
dependency of the cost on the state variables can be very significant.

###### keywords:

SIR model, Optimal Control, Cost-Effectiveness Analysis, Cost Functional,
Infectious disease modeling

††journal: Mathematical Medicine and Biology

##  1 Introduction

Application of the optimal control theory to control problems that arise in
medicine and biology can potentially bring significant benefits, such as
reduction of the cost and side effects of an applied intervention policy.
However, so far, the amount of publications addressing such problems is rather
insufficient and, what is more important, the existing work attracted a lot
less attention to the biologists and medical scientists that they deserved. To
a certain extent, this can be explained by the lack of synergy between the
disciplines and the difficulty of understanding a complicated mathematical
technique. Nevertheless, the authors believe that there is another and,
probably, a more serious reason for the lack of enthusiasm towards the results
of the optimal control theory. The problem is that the results are usually
obtained for a very specific objective functional, and the extents to which
these particular results can be trusted, if this functional is alerted, are
generally unclear. In mathematical terms, it is a problem of robustness and
reliability of the results with respect to perturbations of the objective
functional. While the mathematical theory of robustness for the optimal
control problems is a direction of intensive research, it is far from
completion. Yet, the relevance of the real-life applications urges to explore
this problem using simpler approaches, such as direct computations.

A particular issue that is a very common source of constant controversy in
biomedical applications is a form of the dependency of the cost function on
the considered controls (for example, see discussion and the bibliography in
[22, 4]). The problem is that, in the majority of biomedical applications, the
actual form of the costs and the dependency on these on the controls are
usually unclear and can hardly be defined with a necessary degree of accuracy.
Moreover, the functional dependency can be different under different
circumstances. As a result, for the sake of simplicity, in many publications
in the literature, the cost is postulated to be equal to a sum of the weighted
squares of the considered controls. (Such functionals are typically referred
to as the $L_{2}$-type functionals.) Such a formulation is the most
mathematically convenient, as for such $L_{2}$-type functionals, the optimal
controls can be obtained as explicit functions of the state and adjoint
variables by virtue of the Pontryagin’s maximum principle [21]. Then, the
optimal control problem can be straightforwardly reduced to a two-point
boundary value problem, which is suitable for numerical solving. Thus, the
main advantage of cost functionals with a quadratic term in the control is its
mathematical simplicity. However, in biomedical applications, the
interpretation of such $L_{2}$-type functionals is not always clear and their
use is often difficult to justify. (In engineering, where this type of
objective functional came from, the square of the controls means the energy
spent on the control action [10].)

The problem was properly realized and other forms of dependency were also
considered. For instance, Ledzewicz et al. [13] argued that the $L_{1}$-type
functionals, that is, the cost functionals with a linear dependence on the
control, may capture the cost of an intervention policy more accurately than
the $L_{2}$-type functionals. Di Liddo [4] remarked that, if $u(t)$ is the
fraction of a population $N(t)$ treated by some drug and $p$ is the unitary
cost of the drugs, then a natural way to define the cost $C$ of the treatment
is $C=pu(t)N(t)$. Moreover, as an ultimate way to exclude the controversy of
the intervention cost, Grigorieva et al. [11, 9, 7, 8, 6] suggested to
completely disregard the cost of the intervention under the assumption that it
is ”tolerable” and small compared to the losses that can be inflicted
otherwise. However, one has to take into consideration that the optimal
control with the $L_{1}$-type functionals (and in particular with the cost
depending on the state variables) and the functionals independent of the
controls are generally less amenable to the mathematical analysis, and usually
involves singular and bang-bang controls. Moreover, the linear dependency on
the control might be also not the most realistic option, as some form of the
nonlinearity of the cost is likely to occur, for example, due to the overload
of the health service.

These considerations indicate that further investigation is needed in order to
understand how the objective functional affects the solutions of optimal
control problems, knowledge that can lead to the construction of practical
disease management strategies. Moreover, such a study is essential to verify
the validity and practicality of a large corpus of the results obtained for
$L_{1}$\- and $L_{2}$\- types objective functionals. Another similar related
problem is a form of the functional dependency of the cost on the state
variables and the impact of this on the optimal control. However, despite the
apparent practical relevance of this problem, very little mathematical work
was so far done in this direction: so far, the authors were able to find only
two publications [13, 14] that address this issue.

In this paper, we are addressing this issue considering the problem of the
control of spreading infectious disease as a case study. We consider this
problem for two classical epidemic models, namely a $SIR$ model and a more
complex $SEIRS$ model. For these models, we formulate four optimal control
problems that correspond to four objective functionals that are most frequent
in the literature on disease control, and then analyze the outcomes. We have
to point out that quantitative analysis of the results for different objective
functionals is a nontrivial problem by itself: in this paper, in order to
quantitatively compare the outcomes for different cost functions, we employ
the cost-effectiveness analysis. (The cost-effectiveness analysis was
successfully applied to the epidemiological problem by a number of authors;
for instance, see [1, 2, 12, 16, 17, 20, 23] and bibliography therein.) The
main goal of this study is (i) to investigate how a specific form of the
objective functional impacts the resulting optimal control, and (ii) to
identify if there is a type of objective for which the optimal control
solutions have better performance in terms of cost-efficacy to assess the
advantages or disadvantages of selecting one functional over the others.

The rest of the paper is structured as follows. In the next section, we study
the solutions of an optimal control problem for a classical $SIR$ model with a
single time-dependent control and in section 3 we conduct a numerical study of
the optimal control solutions. In section 4, we carry out a cost-effectiveness
analysis to investigate the cost and health benefits of the control
interventions resulting from the solution of the optimal control problem for
the $SIR$ model. In section 5, we construct a $SEIRS$ model with two time-
dependent controls and formulate an optimal control problem. Section 6
contains a cost-effectiveness analysis based on the incremental cost-
effectiveness ratio (ICER) to measure the performance of the control
interventions for the $SEIRS$ model. Finally, a detailed discussion of the
results is presented in section 7.

##  2 Optimal control for a $SIR$ model with a single control

As a case study, let us consider the classical _SIR_ epidemic model with vital
dynamics given by the following initial value problem:

$\displaystyle\dfrac{dS}{dt}$ $\displaystyle=\mu N-\beta\dfrac{SI}{N}-\mu S,$
(1) $\displaystyle\dfrac{dI}{dt}$ $\displaystyle=\beta\dfrac{SI}{N}-\gamma
I-\mu I,$ $\displaystyle\dfrac{dR}{dt}$ $\displaystyle=\gamma I-\mu R,$

with non-negative initial conditions

$S(0)=S_{0}\geq 0,\quad I(0)=I_{0}\geq 0,\quad R(0)=R_{0}\geq 0.$ (2)

Here, $S(t)$ is the number of susceptibles at time $t$, $I(t)$ the number of
infectious, and $R(t)$ the number of recovered individuals. The total
population $N(t)$ is the sum of the sizes of these three classes, namely,
$N(t)=S(t)+I(t)+R(t)$. In system (1), the parameters $\beta$ and $\gamma$
represent the transmission and the recovery rates, respectively. Individuals
are recruited into the population as susceptibles at a rate $\mu N$. The
natural death of all individuals, whatever their status, occurs also at a per
capita rate $\mu$. Thus, the deaths balance the births, so $N(t)$ is constant
and equal to $N_{0}=N(0)$.

Our main interest here is to identify the optimal healthcare intervention for
curtailing the spread of the infection at the minimal cost. Therefore, we
modify the classical _SIR_ model (1) to include a time-dependent control
function $u(t)$ that represents the treatment of the infectious, and can be
considered as the per capita treatment rate of infected individuals. The _SIR_
model with control becomes

$\displaystyle\dfrac{dS}{dt}$ $\displaystyle=\mu N-\beta\dfrac{SI}{N}-\mu S,$
(3) $\displaystyle\dfrac{dI}{dt}$ $\displaystyle=\beta\dfrac{SI}{N}-\gamma
I-u(t)I-\mu I,$ $\displaystyle\dfrac{dR}{dt}$ $\displaystyle=\gamma
I+u(t)I-\mu R.$

Public health authorities wish to minimize both the cost of infection and the
cost of implementing the control during a time interval $[0,t_{f}]$.
Therefore, we consider the following objective functional:

$J_{i}(u(t))=\int_{0}^{t_{f}}\varphi(I(t))+\sigma_{i}(u(t),I(t))dt.$ (4)

Here, function $\varphi(I)$ represents the cost of the reduction in health and
well-being, that is, the cost of infection. These costs related to pain and
suffering (quality of life) are sometimes referred to as morbidity costs [16].
On the other hand, the function $\sigma_{i}(u,I)$ represents the costs of
healthcare services, that is, the costs of vaccines, treatment,
hospitalizations and other resources.

For simplicity, we assume that the cost of infection is proportional to the
size of the infected class, thus

$\varphi\left(I(t)\right)=A_{1}I(t).$ (5)

The above simplification allows us to focus on the impact of the cost function
$\sigma_{i}(u,I)$ in the form of the optimal control solution. Following the
ideas introduced in [4], we consider the following forms for the cost function
$\sigma_{i}(u,I)$:

$\displaystyle\sigma_{1}(u,I)$ $\displaystyle=$ $\displaystyle
A_{2}u^{2},\quad\text{quadratic state independent cost (QSI)},$ (6)
$\displaystyle\sigma_{2}(u,I)$ $\displaystyle=$ $\displaystyle
A_{2}u^{2}I,\quad\text{quadratic state dependent cost (QSD)}.$ (7)
$\displaystyle\sigma_{3}(u,I)$ $\displaystyle=$ $\displaystyle
A_{2}u,\quad\text{linear state independent cost (LSI)},$ (8)
$\displaystyle\sigma_{4}(u,I)$ $\displaystyle=$ $\displaystyle
A_{2}uI,\quad\text{linear state dependent cost (LSD)}.$ (9)

The coefficients $A_{i}>0$, $i=1,2$ are positive weights on the costs. In
real-life situations, the monetary costs and side effects of a healthcare
intervention or therapy are commonly small compared with the potential losses
that an epidemic can inflict. Therefore, we assume $A_{1}>A_{2}$. In
particular, for later numerical simulations, we will take $A_{1}=100$ and
$A_{2}=10$. We have to remark that, in the opposite case, the monetary costs
of the intervention usually limit in excess the use of a therapy even to the
point of making the policy “not treat any infected people at all” the most
cost-effective strategy [4]. We have to find an optimal control $u^{*}(t)$
such that

$J_{i}(u^{*}(t))=\min_{u\in\mathcal{U}}\left\\{J_{i}(u(t))\right\\}\;\text{subject
to system}\;\eqref{ControlledSIR}$ (10)

where the set of admissible controls is

$\mathcal{U}=\left\\{u:[0,t_{f}]\rightarrow[0,1]\;\text{is Lebesgue
measureble}\right\\}.$ (11)

Let us introduce the vector of states $x=(S,I,R)^{T}$. We call a pair of
states and control $(x,u)$ satisfying both (3) and $u\in\mathcal{U}$ a
feasible pair. Theorem $4.1$ in [5, Chapter III] ensures the existence of an
optimal control $u^{*}(t)$ and corresponding optimal solution
$x^{*}=(S^{*},I^{*},R^{*})^{T}$ for each of our control problems
$\\{QSI,QSD,LSI,LSD\\}$ (the details are given in the Supplementary material
A).

The Pontryagin’s maximum principle states that it is necessary for any optimal
control $u^{*}(t)$ along with the optimal state trajectory
$x^{*}(t)=(S^{*},I^{*},R^{*})^{T}$ to solve the so-called Hamiltonian system.
In other words, this principle converts the optimal control problem (10) into
a problem of minimizing pointwise a Hamiltonian $H$ with respect to the
control $u(t)$. The Hamiltonian is defined in terms of the integrand of the
cost functional $J$ and the right-hand side of the controlled _SIR_ model (3)
as follows:

$\displaystyle H_{i}$
$\displaystyle=A_{1}I+\sigma_{i}(u,I)+\lambda_{S}\left(\mu
N-\beta\dfrac{SI}{N}-\mu S\right)$ (12)
$\displaystyle+\lambda_{I}\left(\beta\dfrac{SI}{N}-\gamma I-u(t)I-\mu
I\right)+\lambda_{R}\left(\gamma I+u(t)I-\mu R\right),$

where $\lambda_{S}$, $\lambda_{I}$ and $\lambda_{R}$ are the associated
adjoint variables.

We apply the Pontryagin’s maximum principle to derive the necessary conditions
that every optimal control satisfies.

######  Theorem 2.1.

Let $u^{*}(t)$ and $x^{*}(t)$ be the optimal control and the corresponding
optimal solutions for the optimal control problem (3), (10). Then there exists
a piecewise differentialble adjoint variable $\lambda(t)$ such that

$H_{i}(t,x^{*},u(t),\lambda(t))\geq H_{i}(t,x^{*}(t),u^{*}(t),\lambda(t))$
(13)

holds at any $t\in[0,t_{f}]$ for all controls $u\in\mathcal{U}$.

Moreover, except at the points of discontinuities of $u^{*}(t)$, the adjoint
variable $\lambda=(\lambda_{S},\lambda_{I},\lambda_{R})^{T}$ satisfies
equations

$\dfrac{d\lambda_{j}}{dt}=-\dfrac{\partial H_{i}}{\partial
x_{j}}(t,x^{*},u^{*},\lambda),\quad j=S,I,R.$ (14)

Therefore,

$\displaystyle\dfrac{d\lambda_{S}}{dt}$
$\displaystyle=\lambda_{S}\left(\beta\dfrac{I^{*}}{N}+\mu\right)-\lambda_{I}\beta\dfrac{I^{*}}{N},$
(15) $\displaystyle\dfrac{d\lambda_{I}}{dt}$
$\displaystyle=-A_{1}-\dfrac{\partial\sigma_{i}}{\partial
I}(u^{*},I^{*})+\lambda_{S}\beta\dfrac{S^{*}}{N}-\lambda_{I}\left(\beta\dfrac{S^{*}}{N}-\gamma-u^{*}-\mu\right)-\lambda_{R}(\gamma+u^{*}),$
$\displaystyle\dfrac{d\lambda_{R}}{dt}$ $\displaystyle=\mu\lambda_{R},$

with transversality conditions

$\lambda_{j}(t_{f})=0,\quad j=S,I,R.$ (16)

The optimality condition is given as

$\dfrac{\partial H_{i}}{\partial u}=\dfrac{\partial\sigma_{i}}{\partial
u}(u,I^{*})-\lambda_{I}I^{*}+\lambda_{R}I^{*}=0.$ (17)

For the quadratic cost functions $\sigma_{i}(u,I)$ ($i=1,2$), using Theorem
2.1, we can straightforwardly obtain the optimality system which consists of
the controlled _SIR_ model, the system (15) for the adjoint variables with
transversality conditions (16), and the characterization of the optimal
control (for the details, see the supplementary material B). The optimality
system allows us to obtain a numerical approximation of the optimal control
for each of this two cases. For the linear cost functions $\sigma_{i}(u,I)$
($i=3,4$), the objective functional is a $L_{1}-$type functional less amenable
to mathematical analysis; therefore, for simplicity, in this case we only
perform a numerical analysis of the optimal controls. The numerical
approximations will allow us to perform the cost-effectiveness analysis.

##  3 Numerical results for the SIR model

We conduct a numerical study of the optimal control problems using BOCOP 2.2.1
software package [3]. In BOCOP, the optimal control problem is approximated by
a finite dimensional optimization problem. This is done by a discretization in
time applied to the state and control variables, as well as the dynamics
equation. These methods are usually less precise than indirect methods based
on the Pontryagin’s maximum principle, but more robust with respect to the
initialization.

(a)  (b)

(c)  (d)

Figure 1: Associated optimal control of the controlled SIR model (3) for the
QSI (a), QSD (b), LSI (c) and LSD (d) cases, respectively.

(a)  (b)

(c)  (d)

Figure 2: Associated optimal states of the controlled SIR model (3) for the
QSI (a), QSD (b), LSI (c) and LSD (d) cases, respectively.

The parameters for the numerical simulations are taken from [19] (Ebola
model): specifically, $\mu=1/(70\times 365)$ (days), $\beta=0.2$
($1/\text{days}$), and $\gamma=0.1$ ($1/\text{days})$. These parameters were
estimated under the assumption that $95\%$ of the population is susceptible
and $5\%$ infectious, therefore, we use initial conditions $N_{0}=100000$,
$S(0)=95000$, $I(0)=5000$, $R(0)=0$. The final time for our simulations is
$t_{f}=100$ days. We consider the optimal control problem (10) denoting the
optimal control functions $u^{*}_{QSI}$, $u^{*}_{QSD}$, $u^{*}_{LSI}$ and
$u^{*}_{LSD}$ for the cases QSI, QSD, LSI and LSD, respectively. The optimal
control profiles are shown in Fig. 2, and the associated optimal states for
each case are shown in Fig. 2.

Fig. 2 indicates that the results can be roughly separated into two groups,
namely, the results for costs QSI and LSI on one hand, and those for QSD and
LSD on another hand. The results for QSD and LSD are very similar. The
similarity between QSI and LSI is not so apparent, but the difference is only
in the rate of the transmission process between the maximum and minimum values
of $u(t)$. (It can be expected that the higher is the value of power $n$ in
the term $u^{n}$ in the objective functional, the slower is the transmission.)
The most remarkable is, however, the difference between these two groups and,
in particular, the visible difference between the costs QSI and QSD (both
these costs are quadratic with respect to $u$), as well as the costs LSI and
LSD (both these costs are linear with respect to $u$). This observation leads
to a conclusion that, at least for the cases where a cost of intervention
policy is low compared with the burden of the disease itself, the optimal
controls for the state-independent intervention costs $A_{2}u^{n}$, for all
$n\in[1,2]$, should be somewhere between the controls depicted in Fig. 2 (a)
and Fig. 2 (c), and that the impact of the actual value of power $n$ in the
cost might be of no principal importance. The same conclusion is even more
apparent for the state-dependent intervention costs of the form $A_{2}u^{n}I$,
where $n\in[1,2]$. That is, the results in the literature that are obtained
for the costs quadratic with respect to control $u(t)$ are valid, at least
qualitatively, for the control problems where costs are linear with respect to
control and for all intermediate costs $u^{n}(t)$ where $n\in[1,2]$. This
implies that results for the objective functions that are quadratic with
respect to controls, which are comparatively easy to deal with mathematically,
are still able to provide a valuable insight into the qualitative behavior of
the optimal controls, and give a correct idea of optimal controls principal
properties. We have to stress that this conclusion is valid only for
comparatively low-cost intervention policies.

At the same time, the results also indicate that an actual form of the
dependency of the intervention costs on the state variables is an issue of
principal importance, and that a variation of this can involve significant
changes in the outcomes; please compare the controls in Fig. 2(a) and Fig.
2(b).

##  4 Cost-effectiveness analysis for the $SIR$ model

Cost-effectiveness analysis is a useful tool that allows to compare the
relative health and economic consequences associated with different strategies
[16]. We compare the cost and health benefits of the control strategies
$\pi_{i}$ defined in terms of the optimal control functions $u_{i}^{*}$ (where
$i=\\{QSI,QSD,LSI,LSD\\}$). In order to do this, we need to compute the cost
and effectiveness for each of these control strategies using some appropriate
health outcome, such as lives saved, the total number of infections averted
and years of life gained.

In this study, to measure the effectiveness of an intervention $\pi$, we
define the cumulative level of averted infections $E(\pi)$ as the difference
between the total numbers of the infected individuals in the absence of
control and of those with the control $\pi$:

$E(\pi)=\int_{0}^{t_{f}}\left(I_{f}(t)-\tilde{I}_{f}(t)\right)dt.$ (18)

Here, $I_{f}(t)$ and $\tilde{I}_{f}(t)$ are, respectively, the numbers of
infected individuals for the no–control scenario and under intervention $\pi$
at time $t$.

To compute the total cost, we can use any of the objective functionals (4).
However, one has to take into consideration the rather obvious fact that if we
use the objective functional $J_{1}$ with cost function $\sigma_{1}(u,I)$
modeling QSI cost, then the corresponding optimal control $u_{QSI}^{*}(t)$ is
the most cost-effective strategy. Therefore, in an attempt to make a fair
comparison of the control strategies for different intervention costs, for
each of the control strategies we compute four costs using all four objective
functionals $J_{i}$, $i=1,2,3,4$. In other words, if $\pi$ is a control
intervention defined by the optimal control function $u(t)$, then the total
cost of applying $\pi$ is defined as

$C_{i}(\pi)=J_{i}(u(t))=\int_{0}^{t_{f}}\varphi(\tilde{I}(t))+\sigma_{i}(u(t),\tilde{I}(t))dt,\quad
i=1,2,3,4.$ (19)

The results of these computations are shown in Table 1.

Strategy | $C_{1}(\pi_{i})$ | $C_{2}(\pi_{i})$ | $C_{3}(\pi_{i})$ | $C_{4}(\pi_{i})$ | $E(\pi_{i})$
---|---|---|---|---|---
$\pi_{QSI}$ | 548971.74 | 603701.05 | 549060.76 | 603701.38 | 800676.1
$\pi_{QSD}$ | 549009.46 | 603691.32 | 549074.40 | 603691.32 | 800676.2
$\pi_{LSI}$ | 549035.59 | 603696.71 | 549036.47 | 603696.71 | 800675.5
$\pi_{LSD}$ | 549009.69 | 603691.32 | 549074.33 | 603691.32 | 800676.2
Table 1: Total costs using each of the cost functions $C_{i}$ ($i=1,2,3,4$)
and effectiveness (infection averted), for the control strategies.

Before we proceed to the cost-effectiveness analysis, we have to point out
that results in Table 1 confirm the above-made conclusion regarding the
qualitative similarity between the optimal controls in strategies QSI and LSI
on one hand, and QSD and LSD on the other hand, and regarding the disparity
between these two couples of scenarios. Indeed, for all four strategies, the
difference between costs $C_{1}(\pi_{i})$ and $C_{3}(\pi_{i})$ in Table 1 is
less than 0.02%, and it is even smaller for the costs $C_{2}(\pi_{i})$ and
$C_{4}(\pi_{i})$. At the same time, for all strategies, the differences
between $C_{1}(\pi_{i})$ and $C_{2}(\pi_{i})$ (and, thus, between
$C_{3}(\pi_{i})$ and $C_{4}(\pi_{i})$) are about 10% .

The cost-effectiveness analysis firstly applies the principle of strong
dominance [16]. According to this principle, if one of the considered
interventions is more effective and less expensive, it is said to be dominant.
If the intervention is more expensive and less effective than an alternative
one, it is said to be dominated. However, for the interventions that are more
effective but at the same time are more costly than the alternatives, the
principle of strong dominance provides no guidance. Please note that, for the
considered model, the effectiveness is practically the same for all four
control strategies (see Table 1). This is due to the fact that the associated
states are almost the same for each of the four considered control functions.
Therefore, the cost-effectiveness analysis is straightforward: the strategy
with the lowest cost in the same time is the most cost-effective strategy.

The corresponding ranking of the strategies is given in Table 2. Please note
that each of the $\pi_{i}$ strategies is the most cost-effective strategy with
respect to its own associated objective functional. This outcome is hardly
unexpected, and implies that the decisive conclusion regarding the cost-
effectiveness of the strategies can be made only on the basis of a specific
cost function, and that if this cost function is used as the objective
functional (as it should be), then the corresponding optimal control and the
optimal solution will be the best-performing strategy.

It can be important to identify the strategy exhibiting the second best
performance. However, due to the fact that the strategies QSI and LSI on one
hand, and QSD and LSD on the other hand, lead to qualitatively similar
outcomes, one can expect that the second best performer in the ranking must
belong to the same couple as the best one. Table 4 confirms this expectation
for three costs of four considered. (Cost $C_{1}$ based of the objective QSI
is the only cost of the four where strategies QSD and LSD exhibited better
cost-efficacy than LSI).

Strategy | $C_{1}-Rank$ | $C_{2}-Rank$ | $C_{3}-Rank$ | $C_{4}-Rank$
---|---|---|---|---
$\pi_{QSI}$ | 1 | 4 | 2 | 4
$\pi_{QSD}$ | 2 | 1 | 4 | 2
$\pi_{LSI}$ | 4 | 3 | 1 | 3
$\pi_{LSD}$ | 3 | 2 | 3 | 1
Table 2: Best-performance ranking for the control strategies using each of the
four possible cost functions. The strategies with the second best performance
are marked in bold.

##  5 Optimal control problem for a $SEIRS$ model with two controls

The above considered $SIR$ model (3) with a single control is one of the
simplest epidemic models in the literature. This simplicity makes the problem
amenable to mathematical analysis; however, at the same time, due to this
simplicity, the influence of different cost functionals on the solution of the
optimal control problem was difficult to detect. Therefore, in this section,
we extend our study to a more complicated case of the compartmental $SEIRS$
epidemic model, which includes two control functions. We believe that analysis
of this model allows to highlight more distinctively the role of a form of the
objective functional on the optimal solutions.

The dynamics of $SEIRS$ model is governed by the following system of
differential equations:

$\displaystyle\dfrac{dS}{dt}$ $\displaystyle=\mu
N-\beta\dfrac{SI}{N}-v(t)S+\theta R-\mu S,$ (20) $\displaystyle\dfrac{dE}{dt}$
$\displaystyle=\beta\dfrac{SI}{N}-\alpha E-\mu E,$
$\displaystyle\dfrac{dI}{dt}$ $\displaystyle=\alpha E-\gamma I-u(t)I-\mu I,$
$\displaystyle\dfrac{dR}{dt}$ $\displaystyle=\gamma I+u(t)I+v(t)S-\theta R-\mu
R,$

where the total population size $N(t)=S(t)+E(t)+I(t)+R(t)$ is assumed to be
constant. The model postulates that the disease has a latent period and,
hence, after infection, individuals pass from the susceptible compartment $S$
into the exposed compartment $E$. The individuals in this compartment are
infected but not yet infectious. The parameter $\alpha$ is the rate at which
the exposed individuals become infectious. The recovered individuals lose the
temporary immunity at a rate $\theta$ and return to the susceptible class $S$.
The rest of parameters are the same as in the $SIR$ model (3).

Compared with the model (3), in this model we incorporate an additional time-
dependent control function $v(t)$ that represents the per capita vaccination
rate of the susceptible individuals. Accordingly, the vaccinated susceptible
individuals enter into the recovered class. We assume that the immunity
acquired through the vaccination is also temporary, and that the vaccinated
susceptible individuals loss immunity at the same rate as the other recovered
individuals. (Such a situation with the loss of immunity occurs if the virus
mutates, as in the case of influenza.) The control problem involves
minimization of the cumulative number of infectious individuals on a given
time interval and of the cost of the intervention policy. The objective
functional is defined as follows:

$J_{i}(v(t),u(t))=\int_{0}^{t_{f}}A_{1}I(t)+\sigma_{i}(v(t),u(t),I(t))dt.$
(21)

Below we consider the following forms for the cost function
$\sigma_{i}(v,u,I)$:

$\displaystyle\sigma_{1}$ $\displaystyle=$ $\displaystyle
A_{2}(v^{2}+u^{2}),\;\text{quadratic state independent cost (QSI)},$ (22)
$\displaystyle\sigma_{2}$ $\displaystyle=$ $\displaystyle
A_{2}(v^{2}S+u^{2}I),\;\text{quadratic state dependent cost (QSD)},$ (23)
$\displaystyle\sigma_{3}$ $\displaystyle=$ $\displaystyle
A_{2}(v+u),\;\text{linear state independent cost (LSI)},$ (24)
$\displaystyle\sigma_{4}$ $\displaystyle=$ $\displaystyle
A_{2}(vS+uI),\;\text{linear state dependent cost (LSD)}.$ (25)

As before, the optimal control problem consists of finding control functions
$v^{*}(t)$ and $u^{*}(t)$, such that

$J_{i}(v^{*},u^{*})=\min_{(v,u)\in\mathcal{U}}\left\\{J_{i}(v(t),u(t))\right\\}\;\text{subject
to system}\;\eqref{ControlledSEIRS}$ (26)

where the set of admissible controls is

$\mathcal{U}=\left\\{(v(t),u(t))|(v(t),u(t))\;\text{measureble},\;0\leq
v(t),u(t)\leq 1,\;t\in[0,t_{f}]\right\\}.$ (27)

The proof of the existence of an optimal control $(v^{*}(t),u^{*}(t))$ and
corresponding optimal solution for problem (26) follows the same argument of
the existence proof for the $SIR$ model with a single control, and is
therefore omitted.

###  5.1 Numerical results for the SEIRS model

(a)  (b)

(c)  (d)

Figure 3: Associated optimal control pair of the controlled SEIRS model (20)
for the QSI (a), QSD (b), LSI (c) and LSD (d) cases, respectively.

(a)  (b)

(c)  (d)

Figure 4: Associated optimal state of the controlled SEIRS model (20) for the
QSI (a), QSD (b), LSI (c) and LSD (d) cases, respectively.

We use the BOCOP 2.2.1 software to obtain numerical approximations of the
solutions to the optimal control problem (26) for each of the objective
functionals. The parameters for the numerical simulations are taken from [24],
specifically, $\mu=0.00003$ ($1/\text{days}$), $\beta=0.25$ ($1/\text{days}$),
$\gamma=0.14$ ($1/\text{days}$), $\alpha=0.33$ ($1/\text{days}$), and
$\theta=0.07$ ($1/\text{days}$). The final time is $t_{f}=100$ days, and the
total population is $N_{0}=100000$.

Figure 4 shows the numerical approximations of the control profiles. Please
note that the optimal controls for the density independent objective
functionals QSI and LSI are similar in qualitative terms. In the same way, the
optimal controls for the density dependent cases QSD and LSD also have very
similar profiles. The associated optimal states are shown in Figure 4.

##  6 Cost-effectiveness analysis for the $SEIRS$ control model

To investigate the performance of the control strategies $\pi_{i}$,
$i=\\{QSD,$ $QSI$, $LSD$, $LSI\\}$, defined in terms of the optimal control
pairs $(v_{i}^{*},u_{i}^{*})$, we perform the cost-effectiveness analysis. As
discussed above, we measure the effectiveness of a control strategy by the
cumulative number of averted infections. To compute the total cost, we use the
objective functionals $J_{i}$ (21). Hence, if $\pi$ is a control intervention
defined by the optimal control pair $(v,u)$, then the total cost of the
intervention $\pi$ is defined as

$C_{i}(\pi)=J_{i}(v,u)=\int_{0}^{t_{f}}A_{1}I+\sigma_{i}(v,u,I)dt,\quad
i=1,2,3,4.$ (28)

To proceed with the comparison of the interventions, we use the model outcomes
to compute the effectiveness and total cost for each of the four considered
control strategies. The results are summarized in Table 3.

Strategy | $C_{1}(\pi_{i})$ | $C_{2}(\pi_{i})$ | $C_{3}(\pi_{i})$ | $C_{4}(\pi_{i})$ | $E(\pi_{i})$
---|---|---|---|---|---
$\pi_{QSI}$ | 1697872.45 | 3836175.45 | 1698013.74 | 5862075.05 | 1267715.57
$\pi_{QSD}$ | 1813853.47 | 2012056.49 | 1813920.59 | 2226189.56 | 1266565.05
$\pi_{LSI}$ | 1697892.72 | 3748946.84 | 1697983.63 | 5155798.23 | 1267715.59
$\pi_{LSD}$ | 1848722.88 | 2032820.69 | 1848757.69 | 2032820.70 | 1266208.63
Table 3: Total costs using each of the cost functions $C_{i}$ ($i=1,2,3,4$)
and effectiveness (infection averted), for the control strategies.

At first glance, from Table 3 one can note that the effectiveness of the
considered strategies does not differ much: the largest difference is of order
of 0.1%. Moreover, $\pi_{QSI}$ and $\pi_{LSI}$ have virtually the same
effectiveness. Furthermore, the costs $C_{1}$ and $C_{3}$ are virtually the
same for all four objective functionals. (The largest difference is smaller
than 0.01%.) At the same time, in contrast to the same data for the $SIR$
model, for the $SEIRS$ model costs $C_{2}$ and $C_{4}$ significantly differ,
in particularly for strategies $\pi_{QSI}$ and $\pi_{LSI}$. (The difference
for $\pi_{QSD}$ is of order of 10%, whereas for $\pi_{LSD}$ costs $C_{2}$ and
$C_{4}$ are virtually the same.) This observation implies that for more
complicated models, and for the objective functional where the control
function appears in the form of a power function $u^{n}$, the cost of the
strategies can depend on the power $n$ more than for a simpler model. However,
the outcomes for objectives QSI and LSI, on one hand, and QSD and LSD, on
another hand, still remain qualitatively similar, and, hence, this dependency
still remains considerably lower than that on the state variables. (We have to
remind that, according to Figs. 4 and 4, the optimal controls themselves
exhibit rather weak dependency on the power $n$, while dependency of the state
variable is significant.)

Using the data in Table 3, we obtain the cost-effectiveness ranking in Table
4. The most important observation from the table is that it confirms the
existence of the stable pairs QSI and LSI, on one hand, and QSD and LSD, on
another hand.

As an illustration of how the cost-effectiveness rankings are computed, we
proceed to show explicitly the calculations for the $C_{4}-Rank$ column of
Table 4. Using the principle of strong dominance, we can rule out intervention
$\pi_{QSI}$, which is dominated by intervention $\pi_{LSI}$ yielding the same
effectiveness at a lower cost.

Strategy | $C_{1}-Rank$ | $C_{2}-Rank$ | $C_{3}-Rank$ | $C_{4}-Rank$
---|---|---|---|---
$\pi_{QSI}$ | 1 | 4 | 2 | 4
$\pi_{QSD}$ | 3 | 1 | 3 | 2
$\pi_{LSI}$ | 2 | 3 | 1 | 3
$\pi_{LSD}$ | 4 | 2 | 4 | 1
Table 4: Best-performance ranking for the control strategies using each of the
four possible cost functions. The strategies with the second best performance
are marked in bold.

However, we have to observe that, for the _SEIRS_ model, the principle of
strong dominance is not sufficient to determine the most cost-effective
intervention, and, therefore, more sophisticated techniques should be applied.
When the more effective strategy is also more costly, the decision maker must
decide if the greater effectiveness justifies the cost of achieving it.
Usually, this can be done by calculating the incremental cost-effectiveness
ratio (ICER) of two (or more) strategies, which is defined as the difference
in the costs divided by the difference in the effectiveness. In mathematical
terms, for two interventions $\pi_{1}$ and $\pi_{2}$,

$ICER(\pi_{1},\pi_{2})=\frac{C(\pi_{2})-C(\pi_{1})}{E(\pi_{2})-E(\pi_{1})},$
(29)

provided that $E(\pi_{1})\neq E(\pi_{2})$. The ICER of each intervention is
found by comparing it to the next most effective intervention. The ICER are
summarized in Table 5).

Strategy | $C_{4}(\pi_{i})$ | $E(\pi_{i})$ | ICER
---|---|---|---
$\pi_{LSI}$ | 5155798.24 | 1267715.59 | $ICER(\pi_{LSI},\pi_{QSD})=2547.48$
$\pi_{QSD}$ | 2226189.56 | 1266565.05 | $ICER(\pi_{QSD},\pi_{LSD})=542.53$
$\pi_{LSD}$ | 2032820.70 | 1266208.63 |
Table 5: Interventions listed in the order of effectiveness together with
their ICER computed by comparing with the next most effective intervention.
($\pi_{QSI}$ is omitted as dominated by $\pi_{LSI}$.)

Alongside to the ICER, we also introduce the average cost-effectiveness ratio
(ACER) as a willingness-to-pay threshold. The ACER is determined with respect
to the “no effect”/“no cost” alternative, and is mathematically given as [18]:

$ACER(s)=\dfrac{C(S)}{E(S)}.$ (30)

As a willingness-to-pay threshold, we use ACER of the less costly
intervention, which in this case is $\pi_{LSD}$ (see Table 5):

$ACER(\pi_{LSD})=\dfrac{2032820.70}{1266208.63}=1.6054.$ (31)

Since $ICER(\pi_{LSI},\pi_{QSD})>\;ACER(\pi_{LSD})$ and
$ICER(\pi_{QSD},\pi_{LSD})>\;ACER(\pi_{LSD})$, the strategy $\pi_{LSD}$ is the
most cost-effective strategy. Similarly, we compute the ACER for $\pi_{QSD}$,

$ACER(\pi_{QSD})=\dfrac{2226189.56}{1266565.05}=1.7576.$ (32)

Since $ICER(\pi_{LSI},\pi_{QSD})>\;ACER(\pi_{LSD})$, the strategy $\pi_{QSD}$
is the second most cost-effective strategy, and thus $\pi_{QSD}$ has Rank #2,
$\pi_{LSI}$ has Rank #3, and $\pi_{QSI}$ has Rank #4.

##  7 Discussion

There are hardly any doubts that the optimal control theory applied to the
problems originated in medicine and biology can potentially bring significant
benefits. At the same time, so far, such applications are rather limited and
generally undervalued. The authors believe that one of the reasons for this is
a controversy with objectives of a control strategy that are formulated in the
form of objective functionals. The problem, as the authors see it, is that the
vast majority of publications in the literature consider objective functionals
that include a sum of weighted squares of the controls. While the functionals
of such a form lead to mathematically convenient problems, their biomedical
interpretation is often uncertain and even dubious. Besides, the objective
functionals where the controls are included as a function $u^{n}$ with $n\neq
2$ (and, in particular, with $n=1$, that appears to be the most natural case
for many biomedical applications) imply the necessity of dealing with a
considerably more complicated mathematical problem and, therefore, are rather
rare in the literature. This controversy leads to a question regarding the
extent to which results obtained from a specific objective functional are
valid and can be trusted. Unfortunately, despite its obvious practical
relevance, so far this issue has attracted rather small attention: so far the
authors were able to find only two publications [13, 14] that address this
issue. Both publications are of the same authors and are dealing with within-
a-patient disease dynamics.

To further address this issue, in particularly concerning infection control at
a level of the entire population, in this paper we explored the impact of a
form of the objective functionals with respect to controls on the structure of
the optimal controls and the associated optimal solutions for two reasonably
simple classical epidemiological models. Specifically, we considered a _SIR_
and a _SEIRS_ models of the spread of an infectious disease in a population.
For these models, we considered the objective functionals where the cost of
control action were assumed to be a linear and a quadratic functions
($L_{2}$\- and $L_{1}$-type functionals). Moreover, for both these types, we
considered the costs that depend of the density of individuals and are
density-independent. To quantitatively compare the impacts of specific forms
of the objective functionals on the optimal controls, in this paper we also
perform cost-effectiveness analysis.

Our numerical simulations show that for both, the $SIR$ model (3) and the
$SEIRS$ model (28), the optimal controls for the both considered density-
independent cases, where the controls are included as a quadratic (the QSI
scenario) and a linear (the LSI scenario) functions, are qualitatively
similar. Likewise, the optimal controls and optimal solutions for both
density-dependent cases, where the controls are included as quadratic (the QSD
scenario) and linear (the LSD scenario) functions, have very similar profiles
as well. At the same time, the optimal controls and the optimal solutions for
the density-independent and the density-dependent cases significantly differ.
This observation allows us to make a very important conclusion that, for the
objective functionals that include the controls in the form of a power
function $u^{n}$, and at least for the cases where the cost of intervention is
notably lower than the potential losses inflicted by the infection (10% in our
cases), the optimal controls and optimal solutions are robust with respect to
the power $n$.

The second conclusion that follows from our results is that the dependency of
the cost of intervention on the state variables (that is, on the densities of
the considered classes) has a remarkably large impact on the optimal control
and the corresponding optimal solution. Hence, for application of the optimal
control theory to real-life practice, a specific form of the dependency of the
costs on the state variables is significantly more important than the
dependency of the cost on the controls themselves. These conclusions indicate
that, for a specific predetermined form of the dependency of the cost of
intervention on the state densities, the most common in the literature results
obtained for the intervention costs that depend on the squares of the controls
are still able to provide valuable insights and can serve as a guideline for
the cases where the real-life dependency of the cost on the controls is
uncertain or is assumed to be a non-quadratic. We have to stress again that
this conclusion is valid for the cases where the cost of intervention is
notably lower than the losses inflicted by the infection. At the same time, we
also have to stress that, while, for a specific predetermined form of the
dependency of the cost on the densities, the optimal controls and optimal
solutions are reasonably robust with respect to a variation of the power $n$,
the actual cost can be notably different for different values of the power
$n$.

## References

  * Agusto and Leite, [2019] Agusto, F. and Leite, M. (2019).  Optimal control and cost-effective analysis of the 2017 meningitis outbreak in nigeria.  Infectious Disease Modelling, 4:161–187.
  * Biswas et al., [2017] Biswas, S., Subramanian, A., ELMojtaba, I. M., Chattopadhyay, J., and Sarkar, R. R. (2017).  Optimal combinations of control strategies and cost-effective analysis for visceral leishmaniasis disease transmission.  PLoS One, 12(2):e0172465.
  * Bonnans et al., [2017] Bonnans, F., Martinon, P., Giorgi, D., Grélard, V., Maindrault, S., Tissot, O., and Liu, J. (2017).  Bocop 2.0. 5—user guide.  URL http://bocop. org, 480.
  * Di Liddo, [2016] Di Liddo, A. (2016).  Optimal control and treatment of infectious diseases. the case of huge treatment costs.  Mathematics, 4(2):21.
  * Fleming and Rishel, [1975] Fleming, W. and Rishel, R. (1975).  Deterministic and stochastic optimal control.
  * Grigorieva et al., [2012] Grigorieva, E., Bondarenko, N., Khailov, E., and Korobeinikov, A. (2012).  Three-dimensional nonlinear optimal control model of wastewater biotreatment.  Neural, Parallel, and Scientific Computations, 20(1):23–36.
  * Grigorieva et al., [2014] Grigorieva, E., Khailov, E., Bondarenko, N., and Korobeinikov, A. (2014).  Modeling and optimal control for antiretroviral therapy.  Journal of Biological Systems, 22(2):199–217.
  * Grigorieva et al., [2013] Grigorieva, E., Khailov, E., and Korobeinikov, A. (2013).  Parametrization of the attainable set for a nonlinear control model of a biochemical process.  Math. Biosci. Eng., 10(4):1067–1094.
  * Grigorieva et al., [2016] Grigorieva, E., Khailov, E., and Korobeinikov, A. (2016).  Optimal control for a sir epidemic model with nonlinear incidence rate.  Math. Model. Nat. Phenom., 11(4):90–105.
  * Grigorieva et al., [2018] Grigorieva, E., Khailov, E., and Korobeinikov, A. (2018).  Optimal control for an seir epidemic model with nonlinear incidence rate.  Studies in Applied Mathematics, 141(3):353–398.
  * [11] Grigorieva, E., Khailov, E., and Korobeinikov, A. (2020a).  Optimal controls of the highly active antiretroviral therapy.  Abstract and Applied Analysis, page 8107106.
  * [12] Grigorieva, E., Khailov, E., and Korobeinikov, A. (2020b).  Optimal quarantine strategies for covid-19 control models.  arXiv preprint arXiv:2004.10614.
  * Ledzewicz et al., [2004] Ledzewicz, U., Brown, T., and Schättler, H. (2004).  Comparison of optimal controls for a model in cancer chemotherapy with l 1-and l 2-type objectives.  Optimization Methods and Software, 19(3-4):339–350.
  * Ledzewicz and Schättler, [2020] Ledzewicz, U. and Schättler, H. (2020).  On the role of the objective in the optimization of compartmental models for biomedical therapies.  Journal of Optimization Theory and Applications, 19(3-4):339–350.
  * Lukes, [1982] Lukes, D. L. (1982).  Differential equations: classical to controlled.  Elsevier.
  * Muennig and Bounthavong, [2016] Muennig, P. and Bounthavong, M. (2016).  Cost-effectiveness analysis in health: a practical approach.  John Wiley & Sons.
  * Okosun et al., [2013] Okosun, K. O., Rachid, O., and Marcus, N. (2013).  Optimal control strategies and cost-effectiveness analysis of a malaria model.  BioSystems, 111(2):83–101.
  * Phillips et al., [2003] Phillips, K. A., Veenstra, D., Bebber, S. V., and Sakowski, J. (2003).  An introduction to cost-effectiveness and cost–benefit analysis of pharmacogenomics.  Pharmacogenomics, 4(3):231–239.
  * Rachah and Torres, [2015] Rachah, A. and Torres, D. F. (2015).  Mathematical modelling, simulation, and optimal control of the 2014 ebola outbreak in west africa.  Discrete Dynamics in Nature and Society, 2015.
  * Rodrigues et al., [2014] Rodrigues, P., Silva, C. J., and Torres, D. F. (2014).  Cost-effectiveness analysis of optimal control measures for tuberculosis.  Bulletin of mathematical biology, 76(10):2627–2645.
  * Saldaña et al., [2019] Saldaña, F., Korobeinikov, A., and Barradas, I. (2019).  Optimal control against the human papillomavirus: Protection versus eradication of the infection.  In Abstract and Applied Analysis, volume 2019. Hindawi.
  * Sharomi and Malik, [2017] Sharomi, O. and Malik, T. (2017).  Optimal control in epidemiology.  Annals of Operations Research, 251(1-2):55–71.
  * Silva and Torres, [2013] Silva, C. J. and Torres, D. F. (2013).  Optimal control for a tuberculosis model with reinfection and post-exposure interventions.  Mathematical Biosciences, 244(2):154–164.
  * Trawicki, [2017] Trawicki, M. B. (2017).  Deterministic seirs epidemic model for modeling vital dynamics, vaccinations, and temporary immunity.  Mathematics, 5(1):7.

Supplemental Materials:

Impact of a cost functional on the optimal control and the cost-effectiveness:
control of a spreading infection as a case study

Fernando Saldaña1, Ariel Camacho2, Andrei Korobeinikov3

1 Instituto de Matemáticas, Campus Juriquilla, 76230, Universidad Nacional
Autónoma de México, Quéretaro, Mexico

2 Facultad de Ciencias, Universidad Autónoma de Baja California, 22860 Baja
California, Mexico

3 School of Mathematics and Information Science, Shaanxi Normal University,
Xián, China

##  Appendix A Existence of the optimal control

Let $x=(S,I,R)^{T}$ and denote the right-hand side of system (3) as the vector
function $F(t,x,u)$. The optimal control problem is

$\min_{u\in\mathcal{U}}J_{i}(u)=\min_{u\in\mathcal{U}}\int_{0}^{t_{f}}L_{i}(t,x,u)dt\quad\text{subject
to system}\;\eqref{ControlledSIR},$ (S1)

where

$L_{i}(t,x,u)=A_{1}I+\sigma_{i}(u,I)$ (S2)

and cost functions $\sigma_{i}(u,I)$, $i=1,2,3,4$, are given by equations
(6)–(9). To prove the existence of an optimal control, we employ Theorem $4.1$
from Fleming and Rishel [5, Chapter III]. This existence theorem states that
the following conditions are sufficient to guarantee the existence of an
optimal control for (S1):

  * (H1)

$F$ is continuous, and there exist positive constants $K_{1}$ and $K_{2}$ such
that

    * (a)

$|F(t,x,u)|\leq K_{1}(1+|x|+|u|)$

    * (b)

$|F(t,\tilde{x},u)-F(t,x,u)|\leq K_{2}|\tilde{x}-x|(1+|u|)$

hold for all $t\in[0,t_{f}]$. Moreover, $F$ can be written as
$F(t,x,u)=\phi(t,x)+\psi(t,x)u$.

  * (H2)

The set of admissible controls $\mathcal{U}$ is closed and convex. Moreover,
there is at least one feasible pair $(x(t),u(t))$ satisfying both (3) and
$u(t)\in\mathcal{U}$.

  * (H3)

$L_{i}(t,x,\cdot)$, $i=1,2,3,4$ is convex on $\mathcal{U}$, and
$L_{i}(t,x,u)\geq g(u)$, where $g$ is continuous and
$|u|^{-1}g(u)\rightarrow+\infty$ as $|u|\rightarrow\infty$, $u\in\mathcal{U}$.

######  Theorem A.1.

Consider the optimal control problem (S1) with control model (3) and cost
functional $J_{i}(u)=\int_{0}^{t_{f}}L_{i}(t,x,u)dt$, $i=1,2,3,4$. Then there
exist an admissible control function $u^{*}(t)\in\mathcal{U}(t_{f})$ such that
$\min J_{i}(u)_{u\in\mathcal{U}}=J_{i}(u^{*})$.

###### Proof.

Since we are considering a constant population, the solution of model (3) are
bounded. In addition, $F\in C^{1}$, and, therefore, (a) and (b) in (H1) are
ensured by suitable bounds on the partial derivatives of $F$ and on
$F(t,0,0)$. Moreover, the state equations are linear with respect to the
controls $u$, and, thus, $F(t,x,u)=\phi(t,x)+\psi(t,x)u$. The existence of a
feasible pair is guaranteed by the Caratheodory theorem [15, pp. 182] for
Cauchy problems. Moreover, for bounded controls on a finite time interval,
$\mathcal{U}$ is clearly closed and convex, and, hence, (H2) is satisfied as
well. The integrand $L_{i}(t,x,\cdot)$ of the cost functional is non-negative
for $i=1,2,3,4$, quadratic for $i=1,2$ and linear for $i=3,4$ with respect to
the controls. Therefore, $L(t,x,\cdot)$ is convex on $\mathcal{U}$.
Furthermore, the set of admissible controls $\mathcal{U}$ is bounded (this
also implies that the term $|u|$ in (H1) can be omitted), and, hence, (H3)
hold vacuously, thus completing the proof. ∎

##  Appendix B Characterization of the optimal control

Here, we characterize the optimal control by direct use of the maximum
principle for the four cost functions (6)–(9). Our aim is to obtain the so-
called optimality system, which consist of the controlled _SIR_ model, the
system (15) for the adjoint variables with transversality conditions (16), and
the characterization of the optimal control. The optimality system allow us to
obtain a numerical approximation of the optimal control $u^{*}(t)$.

###  B.1 Quadratic state-independent cost

Let us consider the quadratic state-independent cost function

$\sigma_{1}(u,I)=A_{2}u^{2}.$

The Hamiltonian is defined by the following equation:

$\displaystyle H_{1}$ $\displaystyle=A_{1}I+A_{2}u^{2}+\lambda_{S}\left(\mu
N-\beta\dfrac{SI}{N}-\mu S\right)$ (S3)
$\displaystyle+\lambda_{I}\left(\beta\dfrac{SI}{N}-\gamma I-u(t)I-\mu
I\right)+\lambda_{R}\left(\gamma I+u(t)I-\mu R\right),$

The system for the adjoint variables is

$\displaystyle\dfrac{d\lambda_{S}}{dt}$
$\displaystyle=\lambda_{S}\left(\beta\dfrac{I^{*}}{N}+\mu\right)-\lambda_{I}\beta\dfrac{I^{*}}{N},$
(S4) $\displaystyle\dfrac{d\lambda_{I}}{dt}$
$\displaystyle=-A_{1}+\lambda_{S}\beta\dfrac{S^{*}}{N}-\lambda_{I}\left(\beta\dfrac{S^{*}}{N}-\gamma-u^{*}-\mu\right)-\lambda_{R}(\gamma+u^{*}),$
$\displaystyle\dfrac{d\lambda_{R}}{dt}$ $\displaystyle=\mu\lambda_{R},$

with transversality conditions
$\lambda_{S}(t_{f})=\lambda_{I}(t_{f})=\lambda_{R}(t_{f})=0$.

The optimality condition for this case is

$\dfrac{\partial H_{1}}{\partial
u}=2A_{2}u-\lambda_{I}I^{*}+\lambda_{R}I^{*}=0.$ (S5)

From the maximum principle we obtain the following characterization for the
optimal control:

$u^{*}(t)=\min\left\\{1,\max\left\\{0,\dfrac{1}{2A_{2}}(\lambda_{I}-\lambda_{R})I^{*}\right\\}\right\\}.$
(S6)

###  B.2 Quadratic state-dependent cost

Here, we consider the quadratic state-dependent cost function

$\sigma_{2}(u,I)=A_{2}u^{2}I.$

The Hamiltonian is given by the following equation:

$\displaystyle H_{2}$ $\displaystyle=A_{1}I+A_{2}u^{2}I+\lambda_{S}\left(\mu
N-\beta\dfrac{SI}{N}-\mu S\right)$ (S7)
$\displaystyle+\lambda_{I}\left(\beta\dfrac{SI}{N}-\gamma I-u(t)I-\mu
I\right)+\lambda_{R}\left(\gamma I+u(t)I-\mu R\right),$

The system for the adjoint variables is

$\displaystyle\dfrac{d\lambda_{S}}{dt}$
$\displaystyle=\lambda_{S}\left(\beta\dfrac{I^{*}}{N}+\mu\right)-\lambda_{I}\beta\dfrac{I^{*}}{N},$
(S8) $\displaystyle\dfrac{d\lambda_{I}}{dt}$
$\displaystyle=-A_{1}-A_{2}u^{2}+\lambda_{S}\beta\dfrac{S^{*}}{N}-\lambda_{I}\left(\beta\dfrac{S^{*}}{N}-\gamma-u^{*}-\mu\right)-\lambda_{R}(\gamma+u^{*}),$
$\displaystyle\dfrac{d\lambda_{R}}{dt}$ $\displaystyle=\mu\lambda_{R},$

with transversality conditions
$\lambda_{S}(t_{f})=\lambda_{I}(t_{f})=\lambda_{R}(t_{f})=0$.

The optimality condition for this case is

$\dfrac{\partial H_{2}}{\partial
u}=2A_{2}uI^{*}-\lambda_{I}I^{*}+\lambda_{R}I^{*}=0.$ (S9)

From the maximum principle we obtain the following characterization for the
optimal control:

$u^{*}(t)=\min\left\\{1,\max\left\\{0,\dfrac{1}{2A_{2}}(\lambda_{I}-\lambda_{R})\right\\}\right\\}.$
(S10)


# First-Principles investigation of the First-order Phase Transition in Kagome
$\mathbf{Na_{2}Ti_{3}Cl_{8}}$

Shuang Liu  Zhi-Bo Yin  Meng-Qiu Long  Yun-Peng Wang<EMAIL_ADDRESS>School of Physics and Electronics, Hunan Key Laboratory for Super-micro
structure and Ultrafast Process, Central South University, 932 South Lushan
Road, Changsha, P. R. China

###### Abstract

The antiferromagnetic kagome material $\mathrm{Na_{2}Ti_{3}Cl_{8}}$ exhibits a
first-order phase transition; the low-temperature phase is characterized by
the trimerizations of Ti ions. In this work, we carry out first-principles
calculations on the crystal and electronic structure of
$\mathrm{Na_{2}Ti_{3}Cl_{8}}$ at the high- and low-temperature phases. The
crystal structures, including the lattice constant and the trimerization of Ti
ions, are well reproduced by taking account a small Hubbard correction of
$1\,\mathrm{eV}$. The calculated total energy landscape reveals a first-order
phase transition with the total energy barrier of $80\,\mathrm{meV}$ per
formula unit. Analysis of the electronic structure indicates a direct metal
bonding among Ti in the trimierized low-temperature phase.

##  I Introduction

The antiferromagnetic kagome $\mathrm{Na_{2}Ti_{3}Cl_{8}}$ material undergoes
a nonpolar-to-polar first-order phase transition as the temperature drops to
around 200 K. Hänni et al. (2017); Kelly et al. (2019) The Ti ions form a
perfect kagome lattice in the high-temperature (HT) phase but trimerize in the
low-temperature (LT) phase. The crystal symmetry is also reduced. The space
group of the HT phase is $R\bar{3}m$ but is $R3m$ in the LT phase. The phase
transition is also characterized by a numerous reduction of the in-plane
lattice constant (from HT $7.4\,\mathrm{\AA}$ to LT $7.0\,\mathrm{\AA}$), and
a substantial reduction ($\sim 9\%$) in the volume from 0.315 to 0.287
$\mathrm{nm}^{3}$ per formula unit. The thermal hysteresis of the phase
transition is about 25 K.Hänni et al. (2017) Both the LT and HT phases are
antiferromagnetic, and a high-spin to low-spin transition occurs along with
the structure transition.amd Gerd Meyer et al. (1995) During the cooling
process, an intermediate ordered phase emerges with only half of Ti ions
participating in the trimerization.Hänni et al. (2017)

The exotic antiferromagnetic order on the kagome lattice is interesting by
itself, and was used as the starting point for understanding the phase
transition in $\mathrm{Na_{2}Ti_{3}Cl_{8}}$.Paul et al. (2020) A trimerized
ground state emerges spontaneously from a Heisenberg antiferromagnet on kagome
lattice.Changlani and Läuchli (2015) The trimerized magnetic order coincides
with the Ti ion trimerization in the $\mathrm{Na_{2}Ti_{3}Cl_{8}}$ material.
In a recent workPaul et al. (2020), the phase transition in
$\mathrm{Na_{2}Ti_{3}Cl_{8}}$ was interpreted as driven by the coupling
between the spin and the lattice degrees of freedom. The rotational symmetry
was proposed to be spontaneously broken by the trimerized magnetic state,
which acts as a precursor driving the structural phase transition.Paul et al.
(2020) Recalling that experimental data unambiguously indicate that the phase
transition in $\mathrm{Na_{2}Ti_{3}Cl_{8}}$ is of first-order. A first-order
phase transition, in general, does not require a spontaneous instability in
the high-temperature phase as a guiding and driving force towards the LT
phase. Instead, a lower free energy of the LT phase is sufficient for
guarantee the phase transition. Theoretical understanding of the first-order
phase transition and the changes of electronic structure after the phase
transition is still lacking.

Figure 1:  The crystal structure of the unit cell of
$\mathrm{Na_{2}Ti_{3}Cl_{8}}$ in (a) the high-temperature phase and the (b)
the low-temperature phase. The (c) perfect and (d) distorted with
trimerizations of the Ti kagome lattice for the high- and low-temperature
phases respectively.

In this work, we provide a first-principles based description of the HT and LT
phases of $\mathrm{Na_{2}Ti_{3}Cl_{8}}$ and the first-order phase transition
between them. The calculated lattice structures, including the lattice
constant and the trimerization of Ti ions, agree with experimental data. The
zero-temperature energy landscape possesses an energy barrier of
$0.08\,\mathrm{eV}$ per formula unit from the HT phase to the LT phases.
Electronic structure analysis reveal the splitting of Ti-$d$ states and a
direct metal-metal bonding in the LT phase.

##  II Computational details

Density functional theory (DFT) calculations are conducted using the Vienna Ab
initio Simulation Package(VASP)Kresse and Hafner (1993, 1994); Kresse and
Furthmüller (1996a, b); Kresse and Joubert (1999). We use the semilocal
generalized gradient approximation (GGA) of Perdew-Burke-Ernzerhof (PBE)Perdew
et al. (1996). The intra-atomic electronic correlation among $d$-orbitals of
Ti ions were taken into account using the DFT+U methodDudarev et al. (1998)
with the Hubbard $U$ parameter set to between 0 and 4 eV. The plane-wave
cutoff energy is set to $400\,\mathrm{eV}$ in all the calculations. The first
Brillouin zone is sampled using a uniform $3\times 3\times 1$ $k$-point mesh.
The crystal structure was optimized until the forces on each atom are smaller
than $0.001\,\mathrm{eV/\AA}$. The solid-state climbing-image nudged elastic
band calculations are performed using the VASP Transition State Tools
(VTST).Henkelman and Johnson (2000); Henkelman et al. (2000); Sheppard et al.
(2012)

##  III Results and Discussions

###  III.1 Crystal structures

The $\mathrm{Na_{2}Ti_{3}Cl_{8}}$ material has a layered structure as depicted
in Fig. 1. The side views show that covalently bonded Ti-Cl layers separated
by Na ions, see Fig. 1(a). For each Ti-Cl layer, one Ti cation is located
within an octahedron formed by six Cl anions. The Cl octahedra are perfect and
the Ti ion is at its center for the HT phase, while the Cl octahedra strongly
distort in the LT phase. The Ti ions of each Ti-Cl layer form a regular kagome
lattice at the HT phase, as shown in Fig. 1(c). However, at the LT phase the
kagome lattice distorts to make Ti ions no longer equal-distanced, see Fig.
1(d).

Figure 2:  The Hubbard U dependence of (a) the lattice constant difference,
(b) the total energy difference between the low-temperature and high-
temperature phases, and (c) the difference between side lengths of Ti-
triangles in the low-temperature phase.

We first study the crystal structures of $\mathrm{Na_{2}Ti_{3}Cl_{8}}$ in the
LT and HT phases. Both the lattice vectors and the atomic coordinates within
the unit cell are allowed to relax. Structural optimizations starts from the
experimental LT and HT structuresHänni et al. (2017). The initial magnetic
configuration is set to the collinear ferromagnetic or the $120^{\circ}$
noncollinear antiferromagnetic configurations. Experimental data suggest an
antiferromagnetic ground state. Applications of the DFT to materials
containing transition-metal elements usually requires an additional and
empirical correction of the electron-electron Coulomb repulsion, which is not
sufficiently accounted for in the local and semilocal density functionals. As
an early transition metal, the Ti element usually requires only a small
Hubbard correction. We therefore perform a series of DFT calculations with the
Hubbard $U$ parameter set to a value between $0.0\,\mathrm{eV}$ and
$4.0\,\mathrm{eV}$.

The calculation results are shown in Fig. 2. The $120^{\circ}$ noncollinear
antiferromagnetic configuration is almost always the ground state; the lattice
constant of the HT phase is always longer than the LT phase; both are in
accord with experiments. We look into three key physical quantities and
compare the calculated values with the experimental data. The first one is the
difference of lattice constant $a$ of the LT and HT phases, defined as $\Delta
a=a(HT)-a(LT)$. We omit the lattice constant $c$ because the trimerization of
Ti ions happens within the $ab$-plane. For the cases with Hubbard $U$ set to
be larger than $2.0\,\mathrm{eV}$, the calculated lattice constant of the HT
phase is only slightly longer ($\sim 0.05\,\mathrm{\AA}$) than the LT phase,
which may be caused by numerical inaccuracy; in fact, the LT structure relaxes
to the HT structure for these cases. Usage of a such larger $U$ shall be
avoided. For the Hubbard $U$ increases from $0.0$ to $2.0\,\mathrm{eV}$, the
calculated $\Delta a$ increases and then decreases. At $U=1\,\mathrm{eV}$, the
calculated $\Delta a$ reaches its maximum of $0.38\,\mathrm{\AA}$
(experimental data: $\sim 0.4\,\mathrm{\AA}$)amd Gerd Meyer et al. (1995);
Hänni et al. (2017); Kelly et al. (2019). The calculated lattice constant $a$
of LT and HT phases matches experimental values in parentheses:
$7.32\,\mathrm{\AA}$ ($7.43\,\mathrm{\AA}$) for the HT phase and
$6.93\,\mathrm{\AA}$ ($6.98\,\mathrm{\AA}$) for the LT phase.amd Gerd Meyer et
al. (1995)

The second key quantity characterizes the trimerization of Ti cations. In the
HT phase, the Ti ions form a perfect kagome lattice, which comprises of two
types of Ti-triangles with opposite orientations. In the LT phase, the side
lengths of one type of Ti-triangles decrease, while the other type of Ti-
triangles expand, resulting in the so-called Ti trimerization. The difference
between side lengths of Ti-triangles (denoted as $\Delta_{\mathrm{Ti-Ti}}$) is
used for characterizing the amplitude of the Ti trimerization. The
experimental value of Ti-Ti distances are $3.071\,\mathrm{\AA}$ and
$3.903\,\mathrm{\AA}$ for the LT phase, thus experimental $\Delta_{\mathrm{Ti-
Ti}}$ is about $0.9\,\mathrm{\AA}$.Hänni et al. (2017) Fig. 2(b) shows that
the calculated $\Delta_{\mathrm{Ti-Ti}}$ is around $1.2\,\mathrm{\AA}$ when
Hubbard $U$ is lower than $1\,\mathrm{eV}$. For $U$=$2\,\mathrm{eV}$ and
larger, the calculated $\Delta_{\mathrm{Ti-Ti}}$ vanishes. Therefore the PBE+U
calculations with $U>2\,\mathrm{eV}$ fails to predict the Ti trimerization in
the LT phase.

The third quantity is the total energy difference: $\Delta
E=E_{\mathrm{HT}}-E_{\mathrm{LT}}$, which shall be positive to guarantee the
LT phase as the ground state. This criteria is satisfied when the Hubbard $U$
is smaller than $1\,\mathrm{eV}$. As the Hubbard $U$ decreases further, the
$\Delta E$ increases to $0.9\,\mathrm{eV}$ per formula unit at
$U=0\,\mathrm{eV}$. The value of $\Delta E$ shall not be too much, otherwise
the HT phase is not possible to appear. This consideration forces us to choose
the Hubbard $U$ to be around $1\,\mathrm{eV}$. Overall, according to the
behaviors of the aforementioned three quantities versus the Hubbard $U$, we
choose $U=1\,\mathrm{eV}$, which reproduces all the experimental observations
fairly well.

###  III.2 The first-order phase transition

Next, we confirm that the phase transition is of first-order by calculating
the total energy landscape. The Wyckoff positions of Ti atoms in the HT phase
are $(0.5,0.5,0)$ but become
$(0.5+\delta_{\mathrm{Ti}},0.5-\delta_{\mathrm{Ti}},0)$ in the LT phase. The
distortion of Ti kagome lattice is characterized by $\delta_{\mathrm{Ti}}$,
which is one of the order parameter for the phase transition:
$\delta_{\mathrm{Ti}}=0$ for the HT phase and $\delta_{\mathrm{Ti}}=0.02$ for
the LT phase, according to experimental data.Hänni et al. (2017) Other
quantities deeply involved in the phase transition is the lattice constants
$a$ and $c$; both exhibit significant changes and hysteresis during the phase
transition.Hänni et al. (2017) The quasi-two-dimensional nature of the
$\mathrm{Na_{2}Ti_{3}Cl_{8}}$ material makes it natural to choose the lattice
constant $a$, but not $c$ as the other order parameter. We construct the total
energy landscape in the order parameter phase spanned by
$\delta_{\mathrm{Ti}}$ and $a$. In practice, we optimize the lattice constant
$c$ and positions of all atoms except Ti, while the lattice constant $a$ and
the fractal positions of Ti ions are fixed.

Figure 3:  (a) The calculated total energy landscape in the
$a$-$\delta_{\mathrm{{Ti}}}$ plane, where $a$ is the lattice constant in unit
of $\mathrm{\AA}$ and $\delta_{\mathrm{{Ti}}}$ is the deviation of Wyckoff
positions (see main text). (b) The calculated minimum energy path between the
low- and high-temperature phases. The total energy of the low-temperature
phase is set to zero as a reference.

The calculated total energy landscape is shown in Fig. 3. The LT phase
corresponds to $\delta_{\mathrm{Ti}}=0.029$ and $a=6.91\,\mathrm{\AA}$ in Fig.
3, while the HT phase corresponds to $\delta_{\mathrm{Ti}}=0.0$ and
$a=7.28\,\mathrm{\AA}$. The total energy landscape in Fig. 3 confirms that
both the LT and the HT phases are located at the local minima of the energy
landscape, hence the transition between them shall be of first-order. One can
read an energy barrier of about $0.8\,\mathrm{eV}$. We also conduct a solid-
state climbing-image nudged elastic bandSheppard et al. (2012) calculation on
the transition from the LT to the HT phase. The total energy relative to that
of the LT phase along the reaction coordinate is plotted in Fig. 3(b). The
transition state, denoted by the star in Fig. 3(b) is higher by
$0.74\,\mathrm{eV}$ than the LT phase, and $0.24\,\mathrm{eV}$ than the HT
phase.

###  III.3 Electronic structures

The Ti ions in $\mathrm{Na_{2}Ti_{3}Cl_{8}}$ have a nominal valence of +2,
thus in the $d^{2}$ configuration. The crystal field of Ti-$d$ electrons due
to Cl ions in the first shell is perfectly octahedral in the HT phase, hence
the Ti-$d$ states split to triply degenerate $t_{\mathrm{2g}}$ orbitals but
only occupied by two $d$ electrons. Therefore the Ti ions exhibit unquenched
orbital moments. The $\mathrm{V^{2+}}$ ions in the $\mathrm{VI_{3}}$ material
exhibit unqueched orbital moments in a similar manner. However the orbital
moments of $\mathrm{Ti^{2+}}$ ions in $\mathrm{Na_{2}Ti_{3}Cl_{8}}$ is
smaller. The spin-only magnetic moment is
$g\sqrt{S(S+1)}\mu_{B}=2.83\,\mu_{B}$ (setting the $g$-factor to 2.0), which
is only slightly larger than experimental value of $2.55\,\mu_{B}$ at
$300\,\mathrm{K}$. This indicates nearly quenched orbital moments for
$\mathrm{Ti}^{2+}$ ions in the HT phase. Our DFT calculations including the
spin-orbit coupling effect confirm that the orbital magnetic moment of
$\mathrm{Ti}^{2+}$ ions is smaller than $0.06\,\mu_{B}$. For the LT phase the
calculated orbital magnetic moments is about $0.01\,\mu_{B}$ per Ti ion
although the Cl-octahedra are strongly distorted. The atomic magnetic moments
thus are dominated by the spin moment, which is $1.51\,\mu_{B}$ and
$0.53\,\mu_{B}$ per Ti ion for the HT and the LT phases according to our DFT
calculations. Although the DFT underestimates the magnetic moments by about
$1.0\,\mu_{B}$ (experimental values are $2.55\,\mu_{B}$ and $1.49\,\mu_{B}$
for HT and LT phases), DFT reproduces the $\sim 1.0\,\mu_{B}$ reduction of
magnetic moments from the HT to the LT phase.

Figure 4:  The calculated electronic density of states of (a) the high-
temperature and (b) the low-temperature phases. The real-space distribution of
the charge density corresponding to the occupied states shown in (a,b) for (c)
the HT and (d) the LT phases; the iso-value of the charge density is set to
$0.06/\mathrm{\AA}^{3}$. Color code for species: green for Cl and blue for Ti;
Na ions are omitted for clarity.

The $\mathrm{Na_{2}Ti_{3}Cl_{8}}$ material is semiconducting and the
calculated band gap is about $0.3\,\mathrm{eV}$ for both the HT and the LT
phases, as shown in Fig. 4(a,b). The application of intra-atomic Coulomb
repulsion makes the triply degenerate $t_{\mathrm{2g}}$ orbitals split. For
the HT phase, the unoccupied DOS has a peak at $0.8\,\mathrm{eV}$ and the
occupied DOS has a peak at $-0.2\,\mathrm{eV}$; their difference of
$1.0\,\mathrm{eV}$ comes from the Hubbard $U$ of the same value. The
additional trimerization distortions in the LT phase split the DOS further, as
indicated by Fig. 4(b). A strong DOS peak is detached and resides at
$-0.7\,\mathrm{eV}$, lower than other valence states. The unoccupied DOS also
shows a clear splitting of about $0.5\,\mathrm{eV}$.

##  IV Conclusions

In summary, we re-investigated the phase transition in kagome material
$\mathrm{Na_{2}Ti_{3}Cl_{8}}$. We confirm the first-order nature of the phase
transition as observed in experiments, and we predict an energy barrier of
0.24 eV per formula unit for this transition. The changes in the lattice
constant and the local magnetic moments are also in accord with experimental
data. The trimerization of Ti ions triggers a direct Ti-Ti bonding.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
(Grant No. 12004439) and computational resources from the High Performance
Computing Center of Central South University.

## References

  * Hänni et al. (2017) N. Hänni, M. Frontzek, J. Hauser, D. Cheptiakov, and K. Krämer, Z. Anorg. Allg. Chem. 643, 2063 (2017).
  * Kelly et al. (2019) Z. A. Kelly, T. T. Tran, and T. M. McQueen, Inorg. Chem. 58, 11941 (2019).
  * amd Gerd Meyer et al. (1995) D. J. H. amd Gerd Meyer, T. Dedecke, and W. Urland, Angew. Chem. Im. Ed. Engl. 34, 71 (1995).
  * Paul et al. (2020) A. Paul, C.-M. Chung, T. Birol, and H. J. Changlani, Phys. Rev. Lett. 124, 167203 (2020).
  * Changlani and Läuchli (2015) H. J. Changlani and A. M. Läuchli, Phys. Rev. B 91, 100407 (2015).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Hafner (1994) G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994).
  * Kresse and Furthmüller (1996a) G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996a).
  * Kresse and Furthmüller (1996b) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996b).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Dudarev et al. (1998) S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev. B 57, 1505 (1998).
  * Henkelman and Johnson (2000) G. Henkelman and D. D. Johnson, J. Chem. Phys. 113, 9978 (2000).
  * Henkelman et al. (2000) G. Henkelman, B. Uberuaga, and D. D. Johnson, J. Chem. Phys. 113, 9901 (2000).
  * Sheppard et al. (2012) D. Sheppard, P. Xiao, W. Chemelewski, D. D. Johnson, and G. Henkelman, J. Chem. Phys. 136, 074103 (2012).