# Semantic Extractor-Paraphraser based Abstractive Summarization

Anubhav Jangra
IIT Patna, India
<EMAIL_ADDRESS>
&Raghav Jain11footnotemark: 1
DTU, India
<EMAIL_ADDRESS>
&Vaibhav Mavi11footnotemark: 1
IIT Delhi, India
<EMAIL_ADDRESS>
Sriparna Saha
IIT Patna, India
<EMAIL_ADDRESS>
&Pushpak Bhattacharyya
IIT Bombay, India
<EMAIL_ADDRESS>
* means equal contribution.

###### Abstract

The anthology of spoken languages today is inundated with textual information,
necessitating the development of automatic summarization models. In this
manuscript, we propose an extractor-paraphraser based abstractive
summarization system that exploits semantic overlap as opposed to its
predecessors that focus more on syntactic information overlap. Our model
outperforms the state-of-the-art baselines in terms of ROUGE, METEOR and word
mover similarity (WMS), establishing the superiority of the proposed system
via extensive ablation experiments. We have also challenged the summarization
capabilities of the state of the art Pointer Generator Network (PGN), and
through thorough experimentation, shown that PGN is more of a paraphraser,
contrary to the prevailing notion of a summarizer; illustrating it’s
incapability to accumulate information across multiple sentences.

##  1 Introduction

Over the past few years, the Internet has become the most convenient and
preferred form of information sharing worldwide. The evolution of technology
has made it possible for anyone to convey their knowledge, opinions and ideals
to the world, resulting in an increasing surge of information hindering users
from accessing desired content. This increasing need to obtain key information
makes the task of summarization paramount. Text is the most widely adopted
form of communication, be it for personal
messaging111https://news.gallup.com/poll/179288/new-era-communication-
americans.aspx or for broadcasting, owing to its ability to convey almost any
concept, its general flexibility to suit everyone’s needs, and its less
storage requirement (opposed to other modes of communication like audio and
video). Text summarization is a problem at the very core of natural language
processing, and has various applications in the spoken languages, including
summarization of conversations, and public speeches.

Some works have been done in the field of reinforcement learning based text
summarization (Dong et al., 2018; Liu et al., 2018), the most prominent
architecture being extractor-abstractor (EXT-ABS) model (Chen and Bansal,
2018). Inspired from this architecture, in this manuscript, we have proposed
an extractor-paraphraser system that uses semantic information overlap as the
underlying guidance strategy. The model is further enhanced to surpass its
limits using reinforcement learning, for which we have proposed a novel
semantic overlap based reward function. Word Mover Similarity (WMS) (Clark et
al., 2019) is utilized to evaluate semantic similarity across generated
sentences and the true ground truth summary sentences.

Figure 1: Proposed model architecture.

We assume that paraphrasing is a relatively simpler task than abstractive
summarization, with the underlying intuition that paraphrasing is a sub-
problem within abstractive summarization. To bolster our hypothesis,
experiments are conducted on the extractor-abstractor (EXT-ABS) model Chen and
Bansal (2018) and the Pointer Generator Network (PGN) See et al. (2017), which
is used as the basic abstraction unit in the former architecture. The results
are rather staggering and reveal that the PGN model also paraphrases input
document sentences, albeit implicitly. The major contributions of the paper
are as follows:

  * •

A novel semantic overlap based reward function is proposed for reinforcement
of extractor-paraphraser model.

  * •

To the best of our knowledge, we are the first ever to discover the fact that
PGN networks are indeed doing an implicit extraction-paraphrasing operation,
revealing the true nature of existing abstractive summarization models.

The rest of this paper is structured as follows: In Section 2 we have
discussed related works of automatic text summarization. In Section 3 we have
described the proposed model, and in Section 4 we have stated the experimental
setup and the datasets used. A thorough discussion and state the results are
provided in Section 5, followed by the conclusion and future work in Section
6.

##  2 Related Work

Automatic text summarization has been extensively researched over more than
three decades, and has shown a lot of progress and promise over the course of
time. Various approaches have been explored to tackle both extractive and
abstractive summarization. Initial research Paice (1990); Kupiec et al. (1995)
focused on extractive summarization due to its easier setup. Various
techniques ranging from integer linear programming Galanis et al. (2012),
graph based approaches Mihalcea and Tarau (2004); Mihalcea (2004), genetic
algorithms Saini et al. (2019a, b), and neural networks Nallapati et al.
(2017); Zhang et al. (2016) have been adopted to solve the extractive
summarization task. The majority of the research in abstractive summarization
revolves around deep learning See et al. (2017); Chopra et al. (2016);
Nallapati et al. (2016). Liu et al. (2018) proposed a generative adversarial
network based model to generate document abstracts. A handful of works however
also use ILP Banerjee et al. (2015) and graph-based Ganesan et al. (2010)
techniques to attempt to solve the problem. A lot of domain specific
summarization techniques have also been explored, like radiology findings
summarization Zhang et al. (2018), across-time summarization Duan and Jatowt
(2019), movie review summarization Zhuang et al. (2006), book summarization
Mihalcea and Ceylan (2007), and customer review based opinion summarization
Pecar (2018). Lately, multi-modal summarization Jangra et al. (2020a, b); Zhu
et al. (2020); Saini et al. (2020) has also gained popularity .

Recently, people have also explored reinforcement learning to tackle the
problem of automatic text summarization in both extractive Dong et al. (2018);
Gao et al. (2019) and abstractive domains Xiao et al. (2020); Chen and Bansal
(2018). Chen and Bansal (2018) have proposed an extractor-abstractor
architecture, separating the relevant data searching part and the paraphrasing
part to individual modules. In this work, we have proposed a system inspired
from Chen and Bansal (2018), stressing on the significance of semantic
information over the traditional syntactic overlap. The literature on text
summarization is rich, and has an abundance of survey papers Yao et al.
(2017); Gambhir and Gupta (2017) to get an in depth overview of the domain.

##  3 Proposed Method

Problem Definition: Given the training data {X, Y} where
$X=\\{d_{1},d_{2},...,d_{N}\\}$ is the set of input documents and
$Y=\\{y_{1},y_{2},...,y_{N}\\}$ is the set of corresponding output summaries,
the task of automatic summarization is defined as the problem of discovering a
function $f:X\mapsto Y$, such that $f(d_{i})=y_{i};\forall
i\in\\{1,2,...,N\\}$.

We have proposed an extractor-paraphraser framework, which is inspired from
the extractor-abstractor (EXT-ABS) framework introduced by Chen and Bansal
(2018). The summarization function $f(.)$ is approximated as the composition,
$f(d_{i})=h(g(d_{i}))$, where the functions $g(.)$ and $h(.)$ are modeled as
the extractor and the paraphraser components of the model, respectively. Given
an input document
$d_{i}=\\{s^{d_{i}}_{1},s^{d_{i}}_{2},...,s^{d_{i}}_{|d_{i}|}\\}$, the
extractor $g(.)$ extracts relevant sentences, acting as the primary noise
filter. These extracted set of sentences are fed to the paraphraser $h(.)$,
which accumulates the information into a concise gist of the extracted
sentences, simulating the natural language generation module in the proposed
system (Fig. 1). The paraphraser in our system is capable of summarizing
multiple extracted sentences to generate one sentence, in contrast with its
predecessor ‘abstractor’ from EXT-ABS model Chen and Bansal (2018) which
rephrases one extracted sentence at a time. The proposed framework consists of
an ‘n-to-one paraphraser’, that compiles $n$ sentences into one sentence,
generating a richer summary. Formally, given a document
$d_{i}=\\{s^{d_{i}}_{1},s^{d_{i}}_{2},...,s^{d_{i}}_{|d_{i}|}\\}$, the
extractor is defined as:

$\begin{split}&g:X\mapsto Z^{|y_{i}|\times n},\\\ &\text{ s.t.
}g(d_{i})=\\{k_{j,l}\,|\,1\leq
k_{j,l}\leq|d_{i}|\\}_{l=1}^{n}{,}_{j=1}^{|y_{i}|}\\\ \end{split}$ (1)

where $k_{p,q}$ represents the index of the $q^{th}$ extracted sentence
corresponding to the $p^{th}$ sentence in the final gold summary. For
modelling the same, we have used an encoder-decoder model, where the encoder
consists of a temporal convolutional model Kim (2014) cascaded with a
bidirectional LSTM network Schuster and Paliwal (1997) and the decoder is a
uni-directional LSTM model Hochreiter and Schmidhuber (1997) (Fig. 1).

The paraphraser function $h(.)$ is defined as follows:

$h(\\{k_{j,l}\\}_{l=1}^{n}{,}_{j=1}^{|y_{i}|},d_{i})=\\{p(\oplus_{l=1}^{n}s^{d_{i}}_{k_{j,l}})\\}_{j=1}^{|y_{i}|}$
(2)

where $\oplus$ represents the concatenation of sentences, $k_{j,l}$ represents
the extractor output and $p$ is written as:

$p(s^{ext}_{j})=s^{y_{i}}_{j}$ (3)

where $s^{ext}_{j}=\oplus_{l=1}^{n}s^{d_{i}}_{k_{j,l}}$ and $s^{y_{i}}_{j}$ is
the $j^{th}$ gold summary sentence. As shown in Fig. 1, a pointer-generator
framework See et al. (2017) is used to model the paraphraser.

###  3.1 Training the Submodules

In order to train the paraphraser beforehand, we need to mimic the output of
an extractor222Note that the extractor’s output can also be used to train the
paraphraser at this step, however, since the extractor’s weights get fine-
tuned using reinforcement learning at a later stage, training paraphraser on
extractor’s non-ideal outputs might lead to unsatisfactory performance of the
paraphraser.. To fulfil this requirement, exemplary-extracted sentences are
generated using gold summaries. Word-mover distance (WMD) Kusner et al. (2015)
is adopted to create these exemplary-extracted sentences, with the motivation
that it captures the semantic overlap between sentences better than its
predecessors. For a training pair $\\{d,y\\}$, the labels are generated as:

$\begin{split}\forall l\in\\{1,2,...,n\\}:\\\ \forall s_{j}\in y:\\\
&k_{j,l}=argmin_{i}\\{WMD(s^{d}_{i},s^{y}_{j})\\};\\\ &\forall s^{d}_{i}\in
d-\\{s^{d}_{k_{j^{\prime},0}}\\}_{j^{\prime}=1}^{j-1}\end{split}$ (4)

where $k_{j,l}$ represents the $l^{th}$ exemplary- extracted sentence
corresponding to the $j^{th}$ gold summary sentence, $s^{y}_{j}$.

The $n$ exemplary-extracted sentences corresponding to each summary sentence
are concatenated and fed into the paraphraser as the input and the summary
sentence is fed as the output. We argue that this should allow the paraphraser
to generate information rich summaries. To facilitate the expected behaviour
during the testing phase, we require the extractor to feed the paraphraser
with input similar to the exemplary-extracted sentences. Hence, we first pre-
train the extractor on these exemplary-extracted sentences as opposed to
random initialization. Cross-entropy loss is used for training the paraphraser
and pre-training the extractor.

###  3.2 Extractor agent

Enhancement using reinforcement learning: Here, the extractor is trained as
part of an actor-critic model Mnih et al. (2016) which takes an action based
on the current state and current value of parameters to maximize a given
reward at each time step, where the action is to extract $n$ sentences,
$\\{k_{t,l}\\}_{l=1}^{n}$, and the state refers to the set of document
sentences, $d_{i}$, and already extracted sentences,
$\\{d_{k_{1,l}},d_{k_{2,l}},...,d_{k_{t-1,l}}\\};\,l\in\\{1,2,...,n\\}$. The
predicted sentences are concatenated and passed to the paraphraser to get an
output sentence. A novel semantic-based reward function using word mover
distance (WMD) is used as the reward function. Since the reward is to be
maximised, WMD needs to be converted to a similarity function, for which, a
generalised version of word mover similarity (WMS), proposed in Clark et al.
(2019), is used. Formally, at a time step $t$, given an action ${j_{t,l}}$,
and a summary sentence, $s_{t}$, the short term reward, $r$ is calculated
as333Note that we obtain the word mover similarity proposed in Clark et al.
(2019) for the case $a=0$, and $b=1$.:

$r=\frac{a+1}{a+e^{b\times WMD(s_{t}^{y},p(s_{t}^{ext}))}}$ (5)

where $s_{t}^{ext}=\oplus_{l=1}^{n}s^{d}_{k_{t,l}}$ , and $\oplus$ represents
the concatenation of sentences, $\\{s^{d},s^{y}\\}$ are the sentences
belonging to a training pair and $a$ and $b$ are the hyper-parameters
introduced444The hyperparameters are set to $a=1$ and $b=0.5$ after extensive
experimentations, and are used throughout this paper..

To avoid redundant phrases and words, tri-gram avoidance through beam search
Paulus et al. (2017) is applied at a sentence level. For a fair comparison,
details regarding the beam search reranking and other nuances of
implementation are kept the same as Chen and Bansal (2018).

Table 1: Evaluation scores for the generated text summary using ROUGE, METEOR and Word Mover Similarity (WMS). ‘Id-ext’ refers to the ideal extractor experiments. The ’-’ denotes unavailability of a score. Models | ROUGE-1 | ROUGE-2 | ROUGE-L | METEOR | WMS
---|---|---|---|---|---
Extractive baselines
$EXT-ABS\,(Ext\,only)$ | 40.17 | 18.11 | 36.41 | 22.81 | -
$O2O\,(Ext\,only)$ | 40.97 | 18.42 | 37.35 | 21.86 | 14.28
$M2O_{n=2}\,(Ext\,only)$ | 40.19 | 17.98 | 36.60 | 21.62 | 14.24
Abstractive baselines
$PGN$ | 39.53 | 17.28 | 36.38 | 18.72 | 13.36
$EXT-ABS\,(w/o\,RL)$ | 38.38 | 16.12 | 36.04 | 19.39 | -
$O2O\,(w/o\,RL)$ | 39.82 | 17.05 | 37.21 | 19.24 | 13.81
$M2O_{n=2}\,(w/o\,RL)$ | 32.82 | 11.29 | 31.08 | 16.13 | 12.92
Reinforced models
$EXT-ABS+RL_{ROUGE}$ | 40.88 | 17.8 | 38.53 | 20.38 | 13.7
$O2O+RL_{ROUGE}$ | 41.32 | 18.16 | 38.89 | 20.52 | 14.56
$EXT-ABS+RL_{WMS}$ | 40.82 | 17.82 | 38.45 | 21.47 | 14.42
$O2O+RL_{WMS}$ | 41.2 | 18.12 | 38.81 | 21.34 | 14.6
$M2O_{n=2}+RL_{ROUGE}$ | 39.71 | 16.7 | 37.32 | 18.25 | 13.52
Ablation experiments (Ideal Extractor)
$Id-ext\,EXT-ABS$ | 49.73 | 26.53 | 47.2 | 24.36 | 19.34
$Id-ext\,O2O$ | 50.04 | 26.31 | 47.34 | 24.61 | 20.14
$Id-ext\,M2O_{n=2}$ | 47.00 | 23.37 | 44.23 | 22.22 | 17.99

##  4 Experiments

###  4.1 Dataset

For all the following experiments we have used the CNN / DailyMail dataset
Nallapati et al. (2016), which contains online news articles, with the bullet
highlights treated as the gold standard summaries. The experiments in this
work are conducted on the non-anonymized version of this dataset. The dataset
consists of 277,226 training, 13,368 validation and 11,490 test article-
summary pairs. An article contains $\sim\\!780$ tokens per document, whereas
the summary consists of $\sim\\!56$ tokens with the average number of
sentences per summary being $\sim\\!3.75$. An article sentence, on average
contains $\sim\\!30$ tokens.

###  4.2 Comparative methods

To highlight the superiority of the proposed semantic-overlap based
methodology over existing syntactic measures, we set the value of $n=1$ for
initial experiments. To compare the complexity of the summarization task with
paraphrasing objectively, the experiments are further extended for $n=2$ as
well 555We have limited our work to $n=2$ since the Two-to-one baselines did
not perform efficaciously. Experiments for $n>2$ would be done in future
works.. We evaluate our model with sufficient baselines 666Statistical
analysis on all the variations of the proposed model has been done.,
including:

Pointer-Generator network [$\mathbf{PGN}$]: An encoder-decoder attention based
framework for abstractive summarization See et al. (2017).

Pre-trained extractor - Extractor-Abstractor framework [$\mathbf{EXT-
ABS\,(Ext\,only)}$]: The pre-trained extractor of the extractor-abstractor
framework proposed by Chen and Bansal (2018) on ROUGE-L based exemplary-
extracted sentences.

Extractor-Abstractor framework without reinforcement [$\mathbf{EXT-
ABS\,(w/o\,RL)}$]: The extractor-abstractor framework proposed by Chen and
Bansal (2018) at a stage before reinforcing the extractor.

Extractor-Abstractor framework [$\mathbf{EXT-}$ $\mathbf{ABS+RL_{X}}$]: The
complete extractor-abstractor framework proposed Chen and Bansal (2018)
including X as the reward function for reinforcement learning (X is either
ROUGE-L or the reward function defined in Eq. 5) along with beam search.

Pre-trained extractor - One-to-one extractor-paraphraser model
[$\mathbf{O2O\,(Ext\,only)}$]: The pre-trained extractor of the proposed
extractor-paraphraser framework (with $n=1$) on WMD based (Eq. 4) exemplary-
extracted sentences.

One-to-one extractor-paraphraser model without reinforcement
[$\mathbf{O2O\,(w/o\,RL)}$]: A particular setting of the proposed methodology
where $n=1$, without any reinforcement or beam search.

One-to-one extractor-paraphraser model [$\mathbf{O2O+}$ $\mathbf{RL_{X}}$]: A
particular setting of the proposed methodology where $n=1$, including
extractor, paraphraser, and X as the reward function for reinforcement
learning (X is either ROUGE-L or the reward function defined in Eq. 5) and
beam search.

Pre-trained extractor - Two-to-one extractor-paraphraser model
[$\mathbf{M2O_{n=2}\,(Ext\,only)}$]: The pre-trained extractor777For fair
comparison of extraction capabilities across all models, we limit the model to
output 4 sentences during evaluation on test dataset. of the proposed
extractor-paraphraser framework (with $n=2$) on WMS based (Eq. 4) exemplary-
extracted sentences.

Two-to-one extractor-paraphraser model without reinforcement
[$\mathbf{M2O_{n=2}\,(w/o\,RL)}$]: A specific case of the proposed many-to-one
paraphrasing where $n=2$. This specific model comprises of only the extractor
and paraphraser modules.

Two-to-one extractor-paraphraser model [$\mathbf{M2O_{n=2}+}$
$\mathbf{RL_{X}}$]: A specific case of the proposed many-to-one paraphrasing
where $n=2$; including X as the reward function for reinforced extractor (X is
either ROUGE-L or the reward function defined in Eq. 5) and beam search.

Ideal extractor-abstractor model [$\mathbf{Id-ext}$ $\mathbf{\,EXT-ABS}$]: To
determine the performance of the abstractor component, the $EXT-ABS$ baseline
is evaluated with the assumption that the extractor is ideal, subsequently
feeding the exemplary-extracted sentences generated for the test data using
ROUGE-L score to the paraphraser directly.

Ideal extractor-paraphraser model [$\mathbf{Id-ext}$ $\mathbf{\,O2O}$]:
Similar to the $Id-ext\,EXT-ABS$ setup, the ’n-to-1’ paraphraser is evaluated
with $n=1$, given that the extractor performs ideally (generating exemplary-
extracted sentences as proposed in Eq. 4).

Ideal extractor-paraphraser model [$\mathbf{Id-ext}$ $\mathbf{\,M2O_{n=2}}$]:
The ’n-to-1’ paraphraser with $n=2$ and the exemplary-extracted sentences
assumed as the extractor (generating exemplary-extracted sentences as proposed
in Eq. 4).

##  5 Results

Results of different baselines and the proposed approach are discussed in this
section. Table 1 illustrates that our proposed techniques perform better than
the rest of the systems. We have used ROUGE-1, ROUGE-2, ROUGE-L Lin (2004),
METEOR Banerjee and Lavie (2005) and word mover similarity (WMS) Clark et al.
(2019) as evaluation metrics. We believe that ROUGE as an evaluation metric is
incapable of judging the quality of an abstractive summary due to its emphasis
on syntactic overlap over semantic overlap Liu et al. (2016); Clark et al.
(2019); Novikova et al. (2017). To overcome this, we have also used WMS as an
evaluation metric.

Figure 2: Average sentence distance (word mover distance) scores for most
similar sentences.

###  5.1 Semantic information overlap

It is noticed that the proposed paraphraser in model $O2O\,$ $(w/o\,RL)$
outperforms the abstractor from $EXT-ABS\,$ $(w/o\,RL)$ in terms of ROUGE
scores, while scoring marginally less in terms of METEOR. The extractor
counterparts $EXT-ABS\,(Ext\,only)$ and $O2O\,(Ext\,only)$ also portray a
similar tendency, with a wider gap in METEOR scores. We also observe that in
the reinforced extractor models, the $O2O+RL_{X}$ setting surpasses the $EXT-
ABS+RL_{X}$ setting in almost all the metrics888Here $X\in\\{ROUGE,WMS\\}$,
which remains same when comparing the two models. A similar trend is also
observed in the ideal extractor experiments, where the $Id-ext\,O2O$ model
beats the $Id-ext\,EXT-ABS$ model in WMS while keeping other evaluation scores
comparable. The above mentioned observations illustrate the true capabilities
of using semantic overlap based exemplary-extracted sentences.

Keeping the main model same as the $EXT-ABS$ framework, and changing the
reward function from ROUGE-L to WMS (Eq. 5), it is observed that the latter
($EXT-ABS\,RL_{WMS}$) attains significantly better METEOR and WMS scores,
while maintaining comparable ROUGE scores. However, the true capabilities of
the WMS reward function (Eq. 5) come into play when we attach it with our
extractor-paraphraser framework; $O2O\,RL_{WMS}$ bests every other models in
terms of WMS, while keeping ROUGE and METEOR comparable with the best attained
scores.

One critical observation is that the reward function introduces a bias in the
evaluation process. It can be clearly observed from the fact that the model
$EXT-ABS+RL_{ROUGE}$ obtains better ROUGE scores while it pales in comparison
to $EXT-ABS+RL_{WMS}$ in terms of WMS. Since we stress that semantic
information overlap is more significant than the syntactic overlap, we believe
that WMS is better suited for the evaluation task as well as a better choice
for the reward function. It is established by the fact that the models using
WMS as the reward function attain comparable ROUGE scores as well (while the
reverse is not true), indicating that incorporating semantic information can
assist in capturing syntactic information as well. Examples of generated
summaries for the $EXT-ABS\,RL_{ROUGE}$ and $O2O\,RL_{WMS}$ models are
illustrated in Fig. 3. An important observation in the generated summaries is
that the former model produces the name ”shao li” which is not present in the
input document or the gold standard summary, whereas this mistake is avoided
by the $O2O\,RL_{WMS}$ model.

###  5.2 Summarization vs paraphrasing

Theoretically the $M2O_{n=2}$ setting should surpass the $O2O$ setting, since
the former has extra input information at the paraphraser stage that the
latter lacks. However, it is observed that this does not happen; in actuality,
the $M2O_{n=2}$ model is outperformed by the $O2O$ model in all aspects (Table
1). After manual scrutiny of exemplary-extracted sentence pairs fed to the
paraphraser and the generated sentences, it is observed that the expected
accumulation of information does not take place. To quantify this observation,
WMD based overlap of information is computed between the $exemplary-
extracted\,sentence$ and the generated sentence. It is discovered that on
average, the more similar sentence has a WMD of 1.775 while the other one
obtains a value of 2.894, illustrating the inability of the paraphraser to
combine information across the two sentences into one.

Hence, we hypothesize that the PGN model intrinsically paraphrases one input
sentence to generate the corresponding sentence in the generated summary,
innately mimicking the extractor paraphraser behaviour. An experiment is
formulated to evince the truth of this proposed hypothesis. For this
experiment, three different collections of documents are used: 1) ground truth
summary, 2) summaries generated by PGN See et al. (2017), and 3) summaries
generated by the $EXT-ABS$ model Chen and Bansal (2018), all corresponding to
the data in the test set. For every summary sentence, its semantic overlap
(using WMD) is computed with every document sentence and distances of the
closest $\alpha$ document sentences are reported in an increasing order. The
results of this experiment can be seen in Fig. 2, and it is noticed that the
average gap between similarities with respect to first and second most similar
sentences is the highest in the case of the PGN model (approx. 1.57), while
the value is around 0.41 for the gold standard summaries. The curve for gold
standard has a steady slope, where as the slope for PGN rises quickly from the
first point to the second point, and then follows a constant growth
afterwards, illustrating a very high similarity with the most similar
sentence, thus endorsing the hypothesis that the PGN intrinsically paraphrases
one input sentence to generate one sentence in final summary.

Figure 3: Example of generated summaries for $EXT-ABS+RL_{ROUGE}$ model,
$O2O+RL_{WMS}$. The text in red denotes the novel information that is not
present in the gold summary, the text in green denotes the information overlap
that is present exclusively in the generated summary and the text in blue
denotes the information covered in all three scenarios.

###  5.3 Error Analysis

As the results illustrate in Table 1, overall the one-to-one setting of the
proposed n-to-one paraphraser performs better than the two-to-one setting,
revealing a major blockade of the existing sequential frameworks - the
incapability of combining multiple sentences into one, deviating from the
ideal notion of ‘abstraction’. We also observe that the extractors outperform
their corresponding $(w/o\,RL)$ counterparts, and are competing with the
complete models as well to some extent, portraying the difficulty of
abstractive summarization over sentence extraction, along with the inability
of current evaluation metrics like ROUGE and METEOR to focus on semantic
information over syntactic subtleties. As elucidated by the ideal extractor
experiments, the paraphraser has great potential for even achieving state-of-
the-art results in the summarization task, provided the extractor works
ideally. Our future work would focus on mollifying the loss introduced at the
extraction step, and obtaining a better harmony in the extractor-paraphraser
network as a whole.

##  6 Conclusion

In this paper, we propose the extractor-paraphraser model, that comprises of
the n-to-one paraphraser and a novel word mover similarity based reward
function. We show that the semantic overlap based techniques surpass the
strong baselines that rely on syntactic information overlap. We also develop
experiments to unveil that the state-of-the-art pointer-generator network
(PGN) indeed paraphrases input sentences intrinsically, and is unable to merge
two sentences into one agglomerate sentence when explicitly conditioned to do
so, changing our perception of sequence-to-sequence abstractive summarization
models. We show that the existing works are still nowhere close to mimicking a
true human summary, portraying the difficulty of true abstraction over its
simplified formulation of extracting and then paraphrasing.

Acknowledgement: Dr. Sriparna Saha would like to acknowledge the support of
Early Career Research Award of Science and Engineering Research Board (SERB)
of Department of Science and Technology, India to carry out this research.

## References

  * Banerjee and Lavie (2005) Satanjeev Banerjee and Alon Lavie. 2005.  Meteor: An automatic metric for mt evaluation with improved correlation with human judgments.  In _Proceedings of the acl workshop on intrinsic and extrinsic evaluation measures for machine translation and/or summarization_ , pages 65–72.
  * Banerjee et al. (2015) Siddhartha Banerjee, Prasenjit Mitra, and Kazunari Sugiyama. 2015.  Multi-document abstractive summarization using ilp based multi-sentence compression.  In _Twenty-Fourth International Joint Conference on Artificial Intelligence_.
  * Chen and Bansal (2018) Yen-Chun Chen and Mohit Bansal. 2018.  Fast abstractive summarization with reinforce-selected sentence rewriting.  _arXiv preprint arXiv:1805.11080_.
  * Chopra et al. (2016) Sumit Chopra, Michael Auli, and Alexander M Rush. 2016.  Abstractive sentence summarization with attentive recurrent neural networks.  In _Proceedings of the 2016 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 93–98.
  * Clark et al. (2019) Elizabeth Clark, Asli Celikyilmaz, and Noah A Smith. 2019.  Sentence mover’s similarity: Automatic evaluation for multi-sentence texts.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 2748–2760.
  * Dong et al. (2018) Yue Dong, Yikang Shen, Eric Crawford, Herke van Hoof, and Jackie Chi Kit Cheung. 2018.  Banditsum: Extractive summarization as a contextual bandit.  _arXiv preprint arXiv:1809.09672_.
  * Duan and Jatowt (2019) Yijun Duan and Adam Jatowt. 2019.  Across-time comparative summarization of news articles.  In _Proceedings of the Twelfth ACM International Conference on Web Search and Data Mining_ , pages 735–743. ACM.
  * Galanis et al. (2012) Dimitrios Galanis, Gerasimos Lampouras, and Ion Androutsopoulos. 2012.  Extractive multi-document summarization with integer linear programming and support vector regression.  In _Proceedings of COLING 2012_ , pages 911–926.
  * Gambhir and Gupta (2017) Mahak Gambhir and Vishal Gupta. 2017.  Recent automatic text summarization techniques: a survey.  _Artificial Intelligence Review_ , 47(1):1–66.
  * Ganesan et al. (2010) Kavita Ganesan, ChengXiang Zhai, and Jiawei Han. 2010.  Opinosis: A graph based approach to abstractive summarization of highly redundant opinions.  In _Proceedings of the 23rd International Conference on Computational Linguistics (Coling 2010)_ , pages 340–348, Beijing, China. Coling 2010 Organizing Committee.
  * Gao et al. (2019) Yang Gao, Christian M Meyer, Mohsen Mesgar, and Iryna Gurevych. 2019.  Reward learning for efficient reinforcement learning in extractive document summarisation.  _arXiv preprint arXiv:1907.12894_.
  * Hochreiter and Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber. 1997.  Long short-term memory.  _Neural computation_ , 9(8):1735–1780.
  * Jangra et al. (2020a) Anubhav Jangra, Adam Jatowt, Mohammad Hasanuzzaman, and Sriparna Saha. 2020a.  Text-image-video summary generation using joint integer linear programming.  In _European Conference on Information Retrieval_ , pages 190–198. Springer.
  * Jangra et al. (2020b) Anubhav Jangra, Sriparna Saha, Adam Jatowt, and Mohammad Hasanuzzaman. 2020b.  Multi-modal summary generation using multi-objective optimization.  In _Proceedings of the 43rd International ACM SIGIR Conference on Research and Development in Information Retrieval_ , page 1745–1748.
  * Kim (2014) Yoon Kim. 2014.  Convolutional neural networks for sentence classification.  _arXiv preprint arXiv:1408.5882_.
  * Kupiec et al. (1995) Julian Kupiec, Jan Pedersen, and Francine Chen. 1995.  A trainable document summarizer.  In _Proceedings of the 18th annual international ACM SIGIR conference on Research and development in information retrieval_ , pages 68–73. ACM.
  * Kusner et al. (2015) Matt Kusner, Yu Sun, Nicholas Kolkin, and Kilian Weinberger. 2015.  From word embeddings to document distances.  In _International conference on machine learning_ , pages 957–966.
  * Lin (2004) Chin-Yew Lin. 2004.  ROUGE: A package for automatic evaluation of summaries.  In _Text Summarization Branches Out_ , pages 74–81, Barcelona, Spain. Association for Computational Linguistics.
  * Liu et al. (2016) Chia-Wei Liu, Ryan Lowe, Iulian V Serban, Michael Noseworthy, Laurent Charlin, and Joelle Pineau. 2016.  How not to evaluate your dialogue system: An empirical study of unsupervised evaluation metrics for dialogue response generation.  _arXiv preprint arXiv:1603.08023_.
  * Liu et al. (2018) Linqing Liu, Yao Lu, Min Yang, Qiang Qu, Jia Zhu, and Hongyan Li. 2018.  Generative adversarial network for abstractive text summarization.  In _Thirty-second AAAI conference on artificial intelligence_.
  * Mihalcea (2004) Rada Mihalcea. 2004.  Graph-based ranking algorithms for sentence extraction, applied to text summarization.  In _Proceedings of the ACL Interactive Poster and Demonstration Sessions_ , pages 170–173.
  * Mihalcea and Ceylan (2007) Rada Mihalcea and Hakan Ceylan. 2007.  Explorations in automatic book summarization.  In _Proceedings of the 2007 joint conference on empirical methods in natural language processing and computational natural language learning (EMNLP-CoNLL)_ , pages 380–389.
  * Mihalcea and Tarau (2004) Rada Mihalcea and Paul Tarau. 2004.  Textrank: Bringing order into text.  In _Proceedings of the 2004 conference on empirical methods in natural language processing_ , pages 404–411.
  * Mnih et al. (2016) Volodymyr Mnih, Adria Puigdomenech Badia, Mehdi Mirza, Alex Graves, Timothy Lillicrap, Tim Harley, David Silver, and Koray Kavukcuoglu. 2016.  Asynchronous methods for deep reinforcement learning.  In _International conference on machine learning_ , pages 1928–1937.
  * Nallapati et al. (2017) Ramesh Nallapati, Feifei Zhai, and Bowen Zhou. 2017.  Summarunner: A recurrent neural network based sequence model for extractive summarization of documents.  In _Thirty-First AAAI Conference on Artificial Intelligence_.
  * Nallapati et al. (2016) Ramesh Nallapati, Bowen Zhou, Caglar Gulcehre, Bing Xiang, et al. 2016.  Abstractive text summarization using sequence-to-sequence rnns and beyond.  _arXiv preprint arXiv:1602.06023_.
  * Novikova et al. (2017) Jekaterina Novikova, Ondřej Dušek, Amanda Cercas Curry, and Verena Rieser. 2017.  Why we need new evaluation metrics for nlg.  _arXiv preprint arXiv:1707.06875_.
  * Paice (1990) Chris D Paice. 1990.  Constructing literature abstracts by computer: techniques and prospects.  _Information Processing & Management_, 26(1):171–186.
  * Paulus et al. (2017) Romain Paulus, Caiming Xiong, and Richard Socher. 2017.  A deep reinforced model for abstractive summarization.  _arXiv preprint arXiv:1705.04304_.
  * Pecar (2018) Samuel Pecar. 2018.  Towards opinion summarization of customer reviews.  In _Proceedings of ACL 2018, Student Research Workshop_ , pages 1–8.
  * Saini et al. (2020) Naveen Saini, Sriparna Saha, Pushpak Bhattacharyya, and Himanshu Tuteja. 2020.  Textual entailment–based figure summarization for biomedical articles.  _ACM Transactions on Multimedia Computing, Communications, and Applications (TOMM)_ , 16(1s):1–24.
  * Saini et al. (2019a) Naveen Saini, Sriparna Saha, Dhiraj Chakraborty, and Pushpak Bhattacharyya. 2019a.  Extractive single document summarization using binary differential evolution: Optimization of different sentence quality measures.  _PloS one_ , 14(11):e0223477.
  * Saini et al. (2019b) Naveen Saini, Sriparna Saha, Anubhav Jangra, and Pushpak Bhattacharyya. 2019b.  Extractive single document summarization using multi-objective optimization: Exploring self-organized differential evolution, grey wolf optimizer and water cycle algorithm.  _Knowledge-Based Systems_ , 164:45–67.
  * Schuster and Paliwal (1997) Mike Schuster and Kuldip K Paliwal. 1997.  Bidirectional recurrent neural networks.  _IEEE transactions on Signal Processing_ , 45(11):2673–2681.
  * See et al. (2017) Abigail See, Peter J. Liu, and Christopher D. Manning. 2017.  Get to the point: Summarization with pointer-generator networks.  _CoRR_ , abs/1704.04368.
  * Xiao et al. (2020) Liqiang Xiao, Lu Wang, Hao He, and Yaohui Jin. 2020.  Copy or rewrite: Hybrid summarization with hierarchical reinforcement learning.  In _AAAI_ , pages 9306–9313.
  * Yao et al. (2017) Jin-ge Yao, Xiaojun Wan, and Jianguo Xiao. 2017.  Recent advances in document summarization.  _Knowledge and Information Systems_ , 53(2):297–336.
  * Zhang et al. (2016) Yong Zhang, Meng Joo Er, Rui Zhao, and Mahardhika Pratama. 2016.  Multiview convolutional neural networks for multidocument extractive summarization.  _IEEE transactions on cybernetics_ , 47(10):3230–3242.
  * Zhang et al. (2018) Yuhao Zhang, Daisy Yi Ding, Tianpei Qian, Christopher D Manning, and Curtis P Langlotz. 2018.  Learning to summarize radiology findings.  _arXiv preprint arXiv:1809.04698_.
  * Zhu et al. (2020) Junnan Zhu, Yu Zhou, Jiajun Zhang, Haoran Li, Chengqing Zong, and Changliang Li. 2020.  Multimodal summarization with guidance of multimodal reference.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pages 9749–9756.
  * Zhuang et al. (2006) Li Zhuang, Feng Jing, and Xiao-Yan Zhu. 2006.  Movie review mining and summarization.  In _Proceedings of the 15th ACM international conference on Information and knowledge management_ , pages 43–50.


11institutetext: Eric Forgoston 22institutetext: Department of Mathematical
Sciences, Montclair State University, 1 Normal Avenue, Montclair, NJ 07043,
USA
22email<EMAIL_ADDRESS>33institutetext: Leah B. Shaw
44institutetext: Department of Applied Science, The College of William & Mary,
P.O. Box 8795, Williamsburg, VA 23187-8795, USA
44email<EMAIL_ADDRESS>55institutetext: Ira B. Schwartz 66institutetext:
Nonlinear Systems Dynamics Section, Plasma Physics Division, Code 6792, U.S.
Naval Research Laboratory, Washington, DC 20375, USA
66email<EMAIL_ADDRESS>
# Inferring unobserved multistrain epidemic sub-populations using
synchronization dynamics

Eric Forgoston  Leah B. Shaw  Ira B. Schwartz

(Received: date / Accepted: date)

###### Abstract

A new method is proposed to infer unobserved epidemic sub-populations by
exploiting the synchronization properties of multistrain epidemic models. A
model for dengue fever is driven by simulated data from secondary infective
populations. Primary infective populations in the driven system synchronize to
the correct values from the driver system. Most hospital cases of dengue are
secondary infections, so this method provides a way to deduce unobserved
primary infection levels. We derive center manifold equations that relate the
driven system to the driver system and thus motivate the use of
synchronization to predict unobserved primary infectives. Synchronization
stability between primary and secondary infections is demonstrated through
numerical measurements of conditional Lyapunov exponents and through time
series simulations.

###### Keywords:

Multistrain disease models; Inferring unobserved populations; Center
manifolds; Synchronization

††journal: Bulletin of Mathematical Biology

##  1 Introduction

Understanding spread of disease requires both observational data and
mathematical modeling of some class or form (Anderson and May, 1991). However,
much of the measured data is quite limited in that it is observed only once
during typically short time intervals and typically is non-stationary. In many
instances of observed disease spread, only case numbers per unit time are
measured. As a result, even for simple, well known diseases, it is difficult
to reconcile the models with the data.

An important component of epidemic modeling that cannot be reconciled with
data are those individuals required to complete the disease path from
susceptible to observed infected, which may consist of asymptomatic,
individuals. Typically, intermittent stages along the disease path in between
susceptible and infected cannot be measured directly. For example, in
childhood diseases, certain susceptible individuals may have come in contact
with someone who is infectious, but the resulting infected individual remains
latent for a period of time. Sub-populations of such intermittent states are
typically not observed directly, and must be inferred (Forgoston and Schwartz,
2013).

Alternative geometric modeling approaches to understanding epidemic spread
have been examined using time series analysis from spatio-temporal case
observations. Tools from nonlinear time series analysis using embedding theory
have been applied to measles data (Schaffer et al., 1993; Blarer and Doebeli,
1999) to examine chaotic-like predictability of case history in the short
term.

Another time series method of analysis for epidemic spread is that of Time-
series Susceptible-Infected-Recovered (TSIR) modeling (Bjornstad et al.,
2002). Local fits of time-series data generate measures of local reproductive
rates of infection. For childhood diseases, the main assumption is that in the
pre-vaccine years, all newborns introduced into the population as susceptible
individuals become infected. However, the model excludes any latency period of
infection, since it only considers models of SIR type. Therefore, the model
does not predict time variations in the latent sub-population.

In general, since epidemic data is limited, detailed modeling of disease
spread is required and thus is widespread. Connections of data with full
models which are higher dimensional are difficult since there may exist many
unobserved sub-populations along the disease path. High dimensional coupled
patch models of cities have limited data time series when compared to the
large simulation model dimension. Such limited data sets imply the need for
accurate lower dimensional models to reduce the number of unknown parameters.
Latency of infection, which introduces a series of exposed classes
approximating mean delay times, is one example that generates high dimensional
models. However, it is known rigorously that the dynamics in higher
dimensional deterministic models often relaxes asymptotically onto lower
dimensional hyper-surfaces, or center manifolds (Schwartz and Smith, 1983;
Shaw et al., 2007). The advantage in doing center manifold reductions is that
if one can only observe certain components of a disease, then it is possible
to explicitly construct a function that relates the unobserved components
(such as latency, or asymptomatic infections) to those explicitly measured or
observed. Such a relationship, in which unobserved sub-populations are related
dynamically to populations of disease that are measured, is closely aligned to
synchronization of coupled systems.

Synchronization has frequently been studied in driver/driven dynamical systems
(also called transmitter/receiver systems), in which information is
transmitted unidirectionally from the driver to the driven system (Boccaletti
et al., 2002). Some studies have focused on controlling parameters in the
driven system until synchronization with the driver is achieved, as a way of
determining unknown parameters in the driver system (Parlitz et al., 1996;
Dedieu and Ogorzalek, 1997; Chen and Lü, 2002; Huang and Lin, 2013). When the
systems are synchronized, observation of the driven system also allows
unobserved variables from the driver system to be observed (Dedieu and
Ogorzalek, 1997). In this paper, we propose a similar approach to determine
unobserved quantities during an epidemic. The driver will be an observed time
series, such as infection prevalence.

Our proposed method to infer unobserved epidemic compartments differs from
previous approaches because we exploit synchronization dynamics rather than
using statistical inference techniques. Previous studies focus more on
integrating over or otherwise accounting for uncertainty in unknown
compartments (Gibson et al., 2004; Lekone and Finkenstädt, 2006). However, it
is frequently seen in epidemic models that the dynamics approach a lower
dimensional center manifold. Examples include susceptible-exposed-infected-
recovered (SEIR) models (Forgoston et al., 2009; Forgoston and Schwartz, 2013)
and multistrain models for dengue fever (Shaw et al., 2007). Such results
suggest that full information about an epidemic system is not needed, and
variables that are more important to the dynamics can be used to deduce other,
unknown variables.

We illustrate with a simple example for which it can be shown analytically
that driving a system of unknown compartments with values from known
compartments will result in convergence to the correct values for the unknown
compartments. Consider a susceptible-infected-recovered (SIR) epidemic in a
population with susceptible, infected, and recovered fractions $s(t)$, $i(t)$,
and $r(t)$, respectively; birth and death rate $\mu$; contact rate $\beta$;
and recovery rate $\sigma$. Suppose that the population dynamics is described
by the system

$\displaystyle\dot{s}=\mu-\beta si-\mu s$ (1a) $\displaystyle\dot{i}=\beta
si-\sigma i-\mu i$ (1b) $\displaystyle\dot{r}=\sigma i-\mu r$ (1c)

If a time series of the infected fraction $i(t)$ is known (e.g., from public
health reporting) but the remaining compartments are unobserved, the
susceptible fraction can be determined by driving a new differential equation
with the known $i(t)$ as follows. Let $s_{d}(t)$ be the driven variable that
we hope will reproduce the correct susceptible fraction, and evolve it
according to

$\dot{s}_{d}=\mu-\beta s_{d}i-\mu s_{d}.$ (2)

The difference between $s$ and $s_{d}$, $\xi=s-s_{d}$, obeys

$\dot{\xi}=-(\beta i+\mu)\xi.$ (3)

Because $\beta i+\mu>0$ for all $t>0$, $\xi$ will approach 0 and $s_{d}$ will
approach $s$ in the long time limit. Thus, given a sufficiently long time
series for $i(t)$, the unknown susceptible levels can be determined for the
latter part of the time series.

The application we will study here is a multistrain model for dengue fever
(Schwartz et al., 2005). Dengue is a mosquito-borne disease with four co-
circulating serotypes. In previous observations, secondary infections led to
more severe illness and constituted the majority of dengue hospital cases
(e.g., Nisalak et al. (2003)), so it is plausible that one may have data about
secondary infections and want to deduce other asymptomatic quantities such as
primary infective levels. It has been shown (Shaw et al., 2007) that the
dynamics of this dengue model lie on a lower dimensional center manifold. In
particular, peaks in primary and secondary infective levels are observed to
coincide (Schwartz et al., 2005). However, the equations for secondary
infectives derived in Shaw et al. (2007) require knowledge of susceptible and
recovered compartments in addition to primary infectives. In this paper, we
will show that primary infective levels can be determined by driving a system
of equations with known secondary infective time series. Such a system has the
practical value of using real measured data to extract unobserved
compartments.

In Section 2, we review the multistrain model for dengue that will be used in
this study. Section 3 presents center manifold analysis showing the
relationship between driven and driver variables. Numerical results are
presented in Section 4 demonstrating the efficacy of the proposed method, and
Section 5 concludes the article.

##  2 The Two-Serotype Model

We begin by describing the compartmental multistrain disease model found in
Schwartz et al. (2005), studying only two co-circulating serotypes for
simplicity, but we remark that we expect the results to hold for arbitrary
numbers of strains. We assume that a given population may be divided into the
following classes that evolve in time:

  1. 1.

Susceptible class $s(t)$ consists of the fraction of the total population that
is susceptible to all serotypes.

  2. 2.

Primary infectious classes $x_{i}(t)$ consist of the fraction of the total
population that is infected with serotype $i$ and is capable of transmitting
serotype $i$ to susceptible individuals.

  3. 3.

Primary recovered classes $r_{i}(t)$ consist of the fraction of the total
population that has recovered from being infected with serotype $i$.

  4. 4.

Secondary infectious classes $x_{ij}(t)$ consist of the fraction of the total
population that is currently infected with serotype $j$, but previously was
infected with serotype $i$, where $i\neq j$.

In this model, susceptible individuals may develop a primary infection from
either of the two serotypes. Upon recovering, the individual is immune to the
strain that caused the primary infection. However, the individual may develop
a secondary infection from the second serotype. The infectiousness of this
secondary infection can be increased through an antibody-dependent enhancement
factor. Upon recovering from the secondary infection, the individual is immune
to both serotypes. A flow diagram for this model is given in Fig. 1.

$\mu$$\boxed{s}$$\beta$$\boxed{x_{1}}$$\boxed{x_{2}}$$\sigma$$\sigma$$\boxed{r_{1}}$$\boxed{r_{2}}$$\beta$$\beta$$\boxed{x_{12}}$$\boxed{x_{21}}$$\sigma$$\phi$$\phi$$\boxed{r_{tot}}$
Figure 1: Flow diagram for a multistrain disease model with two co-circulating
serotypes (Schwartz et al., 2005). Included is the possibility of enhancing
the secondary infectiousness through an antibody-dependent enhancement factor
$\phi$.

The governing equations for the two serotype multistrain disease model are

$\displaystyle\dot{s}=\mu-\beta s(x_{1}+x_{2}+\phi(x_{21}+x_{12})),$ (4a)
$\displaystyle\dot{x}_{1}=\beta s(x_{1}+\phi x_{21})-\sigma x_{1},$ (4b)
$\displaystyle\dot{x}_{2}=\beta s(x_{2}+\phi x_{12})-\sigma x_{2},$ (4c)
$\displaystyle\dot{r}_{1}=\sigma x_{1}-\beta r_{1}(x_{2}+\phi x_{12}),$ (4d)
$\displaystyle\dot{r}_{2}=\sigma x_{2}-\beta r_{2}(x_{1}+\phi x_{21}),$ (4e)
$\displaystyle\dot{x}_{21}=\beta r_{2}(x_{1}+\phi x_{21})-\sigma x_{21},$ (4f)
$\displaystyle\dot{x}_{12}=\beta r_{1}(x_{2}+\phi x_{12})-\sigma x_{12},$ (4g)

where $\mu$ represents a constant birth rate, $\beta$ is the contact rate,
$\phi$ is a constant that determines the antibody-dependent enhancement (ADE),
and $\sigma$ is the rate of recovery, so that $1/\sigma$ is the mean
infectious period. Although the contact rate $\beta$ could be given by a time-
dependent function (e.g. due to seasonal fluctuations in the mosquito vector
population), for simplicity, we assume $\beta$ to be constant. Unlike single
strain models of $SIR$ type with constant contact rate, Eqs. (4a)-(4g) possess
a range of $\phi$ where the endemic equilibrium is unstable.

Rates of infection due to primary infectious individuals have the form $\beta
sx_{i}$, as found in a classical SIR epidemiological model. However, rates of
infection due to secondary infectious individuals are weighted by the ADE
parameter $\phi$, and have the form $\beta\phi sx_{ij}$. If $\phi=1$, then
there is no antibody-dependent enhancement, and the primary and secondary
infectious individuals are equally infectious. If $\phi=2$, then secondary
infectious individuals are twice as infectious as primary infectious
individuals, and so forth. As long as $\phi>1$ the nonlinear terms involving
secondary infectious individuals will contain an antibody-dependent
enhancement factor.

Throughout this article, we use the following parameter values:
$\mu=0.02$(years)-1, $\beta=200$(years)-1, $\phi=3$(years)-1, and
$\sigma=50$(years)-1. These disease parameters are consistent with estimates
previously used in modeling dengue fever and are summarized in Table 1. In
particular, the contact rate $\beta$ corresponds to a reproductive rate of
infection $R_{0}$ of 3.2-4.8, which is consistent with estimates found in
Ferguson et al. (1999) and Nagao and Koelle (2008).

The antibody-dependent enhancement parameter $\phi$ is selected to put the
system in the chaotic regime in which the strains are desynchronized, which is
thought to be the biologically relevant dynamics (c.f. Cummings et al. (2005);
Shaw et al. (2007)). Dengue has four serotypes, but we model only two here for
simplicity. The dynamics are similar for two and four serotypes, although
there are shifts in the locations of bifurcation points and thus of realistic
$\phi$ values (Billings et al., 2007). We anticipate that all the qualitative
results of this paper will hold for more realistic four-serotype models.

Table 1: Model parameters. Parameter |  Value |  Reference
---|---|---
$\mu$, birth rate, years-1 |  $\sim 0.02$ |  Ferguson et al. (1999)
$\beta$, transmission coefficient, years-1 |  $\sim 200$ |  Ferguson et al. (1999)
$\phi$, ADE parameter, years-1 |  $\geq 1$ |  Schwartz et al. (2005)
$\sigma$, recovery rate, years-1 |  $50$ |  Rigau-Perez et al. (1998)
$n$, number of strains |  $2$ |  -

It should be noted that mortality terms have been omitted from Eqs. (4a)-(4g).
In the analysis that follows, it is useful to analytically determine the
endemic steady state. This equilibrium state is not easy to find analytically
when mortality terms are included, but this state is close to the one found
when all mortality occurs after recovery from infection with two serotypes,
and the mortality rate of the other compartments is $\mu_{m}=0$. Previous work
(Shaw et al., 2007) has shown that the dynamics with $\mu_{m}=0$ are
qualitatively similar to the dynamics (and have the same bifurcation
structure) when the mortality rate $\mu_{m}$ is equal to the value of the
birth rate $\mu=0.02$ used in this article. Furthermore, the $\mu_{m}=0$
assumption is physically reasonable since the mortality rate for dengue is low
and the average age at infection is believed to be young (Nisalak et al.,
2003; Shaw et al., 2007)

The governing equations for the two serotype multistrain disease subsystem
that are driven by the secondary infectious individuals of Eqs. (4a)-(4g) are

$\displaystyle\dot{s}_{d}=\mu-\beta s_{d}(x_{d1}+x_{d2}+\phi(x_{21}+x_{12})),$
(5a) $\displaystyle\dot{x}_{1d}=\beta s_{d}(x_{1d}+\phi x_{21})-\sigma
x_{1d},$ (5b) $\displaystyle\dot{x}_{2d}=\beta s_{d}(x_{2d}+\phi
x_{12})-\sigma x_{2d},$ (5c) $\displaystyle\dot{r}_{1d}=\sigma x_{1d}-\beta
r_{1d}(x_{2d}+\phi x_{12}),$ (5d) $\displaystyle\dot{r}_{2d}=\sigma
x_{2d}-\beta r_{2d}(x_{1d}+\phi x_{21}),$ (5e)

where the subscript $d$ signifies that the variable is being driven. As
before, $\mu$ represents a constant birth rate, $\beta$ is the contact rate,
$\phi$ is a constant that determines the antibody-dependent enhancement (ADE),
and $\sigma$ is the rate of recovery, with parameter values shown in Table 1.

Since $r_{1d}$ and $r_{2d}$ are decoupled from Eqs. (4a)-(5c), the center
manifold analysis of the following section will be performed using the simpler
10-dimensional system given by Eqs. (4a)-(5c).

##  3 Center Manifold Analysis

We will reduce the dimension of the system given by Eqs. (4a)-(5c) using the
center manifold of the system. The analysis begins by determining the endemic
equilibrium state of the system. It is given as

$(s_{0},x_{i,0},r_{i,0},x_{ij,0},s_{d,0},x_{id,0})=\left(\frac{\sigma}{\beta(1+\phi)},\frac{\mu}{2\sigma},\frac{\sigma}{\beta(1+\phi)},\frac{\mu}{2\sigma},\frac{\sigma}{\beta(1+\phi)},\frac{\mu}{2\sigma}\right)$
(6)

for all $i,j$.

A general nonlinear system may be transformed so that the system’s linear part
has a block diagonal form consisting of three matrix blocks. The first matrix
block will possess eigenvalues with positive real part; the second matrix
block will possess eigenvalues with negative real part; and the third matrix
block will possess eigenvalues with zero real part. These three matrix blocks
are respectively associated with the unstable eigenspace, the stable
eigenspace, and the center eigenspace. If there are no eigenvalues with
positive real part, then the orbits will rapidly decay to the center
eigenspace.

Equations (4a)-(5c) cannot be written in a block diagonal form with one matrix
block possessing eigenvalues with negative real part and the other matrix
block possessing eigenvalues with zero real part. Even though it is possible
to construct a center manifold from a system not in separated block form
(Chicone and Latushkin, 1997), it is much easier to apply the center manifold
theory to a system with separated stable and center directions. Therefore, we
transform the original system given by Eqs. (4a)-(5c) to a new system of
equations that will have the eigenvalue structure that is needed to apply
center manifold theory. The theory allows one to find an invariant center
manifold that passes through a fixed point and to which one can restrict the
new transformed system.

###  3.1 Transformation of the Two-Serotype Model

To ease the analysis, we define a new set of variables, $\bar{s}$,
$\bar{x}_{i}$, $\bar{r}_{i}$, $\bar{x}_{ij}$, $\bar{s}_{d}$, and
$\bar{x}_{id}$ for all $i,j$ as $\bar{s}(t)=s(t)-s_{0}$,
$\bar{x}_{i}(t)=x_{i}(t)-x_{i,0}$, $\bar{r}_{i}(t)=r(t)-r_{i,0}$,
$\bar{x}_{ij}(t)=x_{ij}(t)-x_{ij,0}$, $\bar{s}_{d}(t)=s_{d}(t)-s_{d,0}$,
$\bar{x}_{id}(t)=x_{id}(t)-x_{id,0}$. These new variables are substituted into
Eqs. (4a)-(5c).

Then, treating $\mu$ as a small parameter, we rescale time by letting
$t=\mu\tau$. We may then introduce the following rescaled parameters:
$\beta=\beta_{0}/\mu$ and $\sigma=\sigma_{0}/\mu$, where $\beta_{0}$ and
$\sigma_{0}$ are $\mathcal{O}(1)$. The inclusion of the parameter $\mu$ as a
new state variable means that the terms in our rescaled system which contain
$\mu$ are now nonlinear terms. Furthermore, the system is augmented with the
auxiliary equation $\frac{d\mu}{d\tau}=0$. The addition of this auxiliary
equation contributes an extra simple zero eigenvalue to the system and adds
one new center direction that has trivial dynamics. The shifted and rescaled,
augmented system of equations now have the endemic fixed point located at the
origin and are provided in Appendix A.

The Jacobian of these shifted and rescaled equations (Eqs. (21a)-(21k)) is
computed to zeroth-order in $\mu$ and is evaluated at the origin. Ignoring the
$\mu$ components, the Jacobian has only eight linearly independent
eigenvectors. Therefore, the Jacobian is not diagonalizable. However, it is
possible to transform Eqs. (21a)-(21j) to a block diagonal form with a
separated eigenvalue structure. As mentioned previously, this block structure
makes the center manifold analysis easier. We use a transformation matrix,
${\bf P}$, consisting of the eight linearly independent eigenvectors of the
Jacobian along with two other vectors chosen to be linearly independent. There
are many choices for these ninth and tenth vectors; our choice is predicated
on keeping the vectors as simple as possible. This transformation matrix is
given as follows:

${\bf P}=\left[\begin{array}[]{cccccccccc}0&0&0&0&1&0&0&0&0&0\\\
-\phi&0&0&0&0&1&0&0&0&0\\\ 0&-\phi&0&0&0&0&1&0&0&0\\\
\phi&0&0&0&0&0&0&1&0&0\\\ 0&\phi&0&0&0&0&0&0&1&0\\\ 1&0&0&0&0&1&0&0&0&0\\\
0&1&0&0&0&0&1&0&0&0\\\ 0&0&1&0&0&0&0&0&0&1\\\ -\phi&0&\phi&-1&0&1&0&0&0&0\\\
0&-\phi&0&1&0&0&1&0&0&0\\\ \end{array}\right].$ (7)

Using the fact that
$(\bar{s},\bar{x}_{1},\bar{x}_{2},\bar{r}_{1},\bar{r}_{2},\bar{x}_{21},\bar{x}_{12},\bar{s}_{d},\bar{x}_{1d},\bar{x}_{2d})^{T}={\bf
P}\cdot{\bf W}^{T}$, where ${\bf
W}=(W_{1},W_{2},W_{3},W_{4},W_{5},W_{6},W_{7},W_{8},W_{9},W_{10})$, then the
transformation matrix leads to the definition of new variables, $W_{i}$,
$i=1\ldots 10$ that can be found in Appendix B. The application of the
transformation matrix to Eqs. (21a)-(21j) leads to a set of transformed
evolution equations that are found in Appendix C.

###  3.2 Application of the Center Manifold Theory

The Jacobian of Eqs. (23a)-(23k) to zeroth-order in $\mu$ and evaluated at the
origin is

$\left[\begin{array}[]{cccc|cccccc}-\sigma_{0}&0&0&0&0&0&0&0&0&0\\\
0&-\sigma_{0}&0&0&0&0&0&0&0&0\\\
0&0&-\frac{\phi\sigma_{0}}{1+\phi}&0&0&0&0&0&0&0\\\
0&0&0&-\frac{\phi\sigma_{0}}{1+\phi}&0&0&0&0&0&0\\\ \hline\cr
0&0&0&0&0&-\sigma_{0}&-\sigma_{0}&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&\sigma_{0}&-\sigma_{0}&0&0&0\\\
0&0&0&0&0&-\sigma_{0}&\sigma_{0}&0&0&0\\\
0&0&0&0&0&-\sigma_{0}&-\sigma_{0}&0&0&0\\\ \end{array}\right],$ (8)

which shows that Eqs. (23a)-(23k) may be rewritten in the form

$\displaystyle\frac{d{\bf x}}{d\tau}={\bf A}{\bf x}+{\bf f}({\bf x},{\bf
y},\mu),$ (9) $\displaystyle\frac{d{\bf y}}{d\tau}={\bf B}{\bf y}+{\bf g}({\bf
x},{\bf y},\mu),$ (10) $\displaystyle\frac{d\mu}{d\tau}=0,$ (11)

where ${\bf x}=(W_{1},W_{2},W_{3},W_{4})$, ${\bf
y}=(W_{5},W_{6},W_{7},W_{8},W_{9},W_{10})$, ${\bf A}$ is a constant matrix
with eigenvalues that have negative real parts, ${\bf B}$ is a constant matrix
with eigenvalues that have zero real parts, and ${\bf f}$ and ${\bf g}$ are
nonlinear functions in ${\bf x}$, ${\bf y}$ and $\mu$. In particular, ${\bf
A}$ is the upper left block matrix of Eq. (8), while ${\bf B}$ is the lower
right block matrix of Eq. (8).

Therefore, this new system of equations, which is an exact transformation of
Eqs. (4a)-(5c), will rapidly collapse onto a lower-dimensional manifold given
by center manifold theory (Carr, 1981; Chicone and Latushkin, 1997; Duan et
al., 2003). Furthermore, the lower dimensional center manifold is given by

$W_{i}=h_{i}(\bm{y},\mu)$ (12)

where $h_{i}$, $i=1\ldots 4$ are unknown functions.

Substitution of the center manifold functions $W_{i}=h_{i}$ given by Eq. (12)
into the transformed evolution equations given in Appendix C leads to the
center manifold condition given in Appendix D.

In general, it is not possible to solve the center manifold condition for the
four unknown functions $h_{i}(\bm{y},\mu)$, $i=1\ldots 4$. Therefore, a Taylor
series expansion of $h_{i}(\bm{y},\mu)$, $i=1\ldots 4$ in $\bm{y}$ and $\mu$
is substituted into the four equations that comprise the center manifold
condition (Eqs. (24a)-(24d)). The unknown coefficients are determined by
equating terms of the same order, and the center manifold equations are found
to be

$\displaystyle W_{1}=$
$\displaystyle\frac{-\beta_{0}}{\sigma_{0}}W_{5}W_{6}-2\frac{\beta_{0}}{\sigma_{0}}W_{6}W_{7}+\frac{\beta_{0}}{\sigma_{0}}W_{6}W_{9}+\mathcal{O}(\epsilon^{3}),$
(13a) $\displaystyle W_{2}=$
$\displaystyle\frac{-\beta_{0}}{\sigma_{0}}W_{5}W_{7}-2\frac{\beta_{0}}{\sigma_{0}}W_{6}W_{7}+\frac{\beta_{0}}{\sigma_{0}}W_{7}W_{8}+\mathcal{O}(\epsilon^{3}),$
(13b) $\displaystyle W_{3}=$
$\displaystyle-\frac{\beta_{0}(1+\phi)^{2}}{\phi^{2}\sigma_{0}}W_{5}W_{6}-\frac{\beta_{0}(1+\phi)^{2}}{\phi^{2}\sigma_{0}}W_{5}W_{7}+\frac{\beta_{0}(1+\phi)^{2}}{\phi^{2}\sigma_{0}}W_{6}W_{10}$
$\displaystyle+\frac{\beta_{0}(1+\phi)^{2}}{\phi^{2}\sigma_{0}}W_{7}W_{10}+\mathcal{O}(\epsilon^{3}),$
(13c) $\displaystyle W_{4}=$
$\displaystyle-\frac{\beta_{0}(1+\phi)^{2}}{\phi\sigma_{0}}W_{5}W_{7}+\frac{\beta_{0}(1+\phi)^{2}}{\phi\sigma_{0}}W_{7}W_{10}+\mathcal{O}(\epsilon^{3}),$
(13d)

where $\epsilon=|(W_{5},W_{6},W_{7},W_{8},W_{9},W_{10},\mu)|$ so that
$\epsilon$ provides a count of the number of $W_{5}$, $W_{6}$, $W_{7}$,
$W_{8}$, $W_{9}$, $W_{10}$, and $\mu$ factors in any one term. It is worth
noting that the center manifold equations may equivalently be found using a
normal form coordinate transform. Details of the method can be found in
Roberts (2008) and application to epidemic models can be found in Forgoston et
al. (2009) and Forgoston and Schwartz (2013)

Using the relations between $W_{i}$ variables and the “barred” variables given
by Eq. (22), Eqs. (13a)-(13d) can be equivalently written as follows:

$\displaystyle\sigma_{0}(\bar{x}_{1}-\bar{x}_{21})=\beta_{0}(\bar{x}_{1}+\phi\bar{x}_{21})\left[\bar{s}-\bar{r}_{2}+\frac{3\phi}{(1+\phi)}\bar{x}_{12}+\frac{2-\phi}{(1+\phi)}\bar{x}_{2}\right],$
(14a)
$\displaystyle\sigma_{0}(\bar{x}_{2}-\bar{x}_{12})=\beta_{0}(\bar{x}_{2}+\phi\bar{x}_{12})\left[\bar{s}-\bar{r}_{1}+\frac{3\phi}{(1+\phi)}\bar{x}_{21}+\frac{2-\phi}{(1+\phi)}\bar{x}_{1}\right],$
(14b)
$\displaystyle\sigma_{0}\left[\bar{x}_{d1}+\bar{x}_{d2}-\bar{x}_{1}-\bar{x}_{2}\right]=\frac{\beta_{0}(1+\phi)^{2}}{\phi^{2}}\left[\bar{x}_{1}+\bar{x}_{2}-\bar{x}_{d1}-\bar{x}_{d2}-\phi\bar{s}+\phi\bar{s}_{d}\right]\times$
$\displaystyle\hskip
113.81102pt\left[\frac{\bar{x}_{1}+\bar{x}_{2}+\phi\bar{x}_{12}+\phi\bar{x}_{21}}{(1+\phi)}\right],$
(14c)
$\displaystyle\sigma_{0}\left[\bar{x}_{d2}-\bar{x}_{2}\right]=\frac{\beta_{0}(1+\phi)^{2}}{\phi^{2}}\left[\bar{x}_{1}+\bar{x}_{2}-\bar{x}_{d1}-\bar{x}_{d2}-\phi\bar{s}+\phi\bar{s}_{d}\right]\times\left[\frac{\bar{x}_{2}+\phi\bar{x}_{12}}{(1+\phi)}\right].$
(14d)

It is possible to simplify Eq. (14c) by substituting Eq. (14d) into Eq. (14c).
The resulting simplified Eq. (14c) is

$\sigma_{0}\left[\bar{x}_{d1}-\bar{x}_{1}\right]=\frac{\beta_{0}(1+\phi)^{2}}{\phi^{2}}\left[\bar{x}_{1}+\bar{x}_{2}-\bar{x}_{d1}-\bar{x}_{d2}-\phi\bar{s}+\phi\bar{s}_{d}\right]\times\left[\frac{\bar{x}_{1}+\phi\bar{x}_{21}}{(1+\phi)}\right].$
(15)

Solving Eqs. (14a)-(14b) for $\bar{x}_{1}$, $\bar{x}_{2}$ leads to the
following approximation for the invariant manifold onto which the driver
system collapses:

$\displaystyle\bar{x}_{1}=\frac{\sigma_{0}\bar{x}_{21}+\beta_{0}\phi\bar{x}_{21}\left[\bar{s}-\bar{r}_{2}+\frac{2-\phi}{1+\phi}\bar{x}_{2}+\frac{3\phi}{1+\phi}\bar{x}_{12}\right]}{\sigma_{0}-\beta_{0}\left[\bar{s}-\bar{r}_{2}+\frac{2-\phi}{1+\phi}\bar{x}_{2}+\frac{3\phi}{1+\phi}\bar{x}_{12}\right]}$
(16a)
$\displaystyle\bar{x}_{2}=\frac{\sigma_{0}\bar{x}_{12}+\beta_{0}\phi\bar{x}_{12}\left[\bar{s}-\bar{r}_{1}+\frac{2-\phi}{1+\phi}\bar{x}_{1}+\frac{3\phi}{1+\phi}\bar{x}_{21}\right]}{\sigma_{0}-\beta_{0}\left[\bar{s}-\bar{r}_{1}+\frac{2-\phi}{1+\phi}\bar{x}_{1}+\frac{3\phi}{1+\phi}\bar{x}_{21}\right]}.$
(16b)

This invariant manifold was found previously in Shaw et al. (2007) and
represents the close relationship of primary infectives and secondary
infectives.

Similarly, solving the simplified Eq. (14c) given by Eq. (15) and Eqs. (14d)
simultaneously for $\bar{x}_{d1}$ and $\bar{x}_{d2}$ leads to an approximate
invariant manifold for the driven primary infectives expressed almost entirely
in terms of driver system variables:

$\displaystyle\bar{x}_{d1}=\frac{\left(\sigma_{0}\bar{x}_{1}\left(1+\phi\right)+\Gamma\left[\bar{x}_{1}\left(\bar{x}_{1}+\bar{x}_{2}+\phi\left(\bar{x}_{12}+\bar{x}_{21}\right)\right)+\phi\left(\bar{s}_{d}-\bar{s}\right)\left(\bar{x}_{1}+\phi\bar{x}_{21}\right)\right]\right)}{\sigma_{0}\left(1+\phi\right)+\Gamma\left(\bar{x}_{1}+\bar{x}_{2}+\phi\left(\bar{x}_{12}+\bar{x}_{21}\right)\right)}$
(17a)
$\displaystyle\bar{x}_{d2}=\frac{\left(\sigma_{0}\bar{x}_{2}\left(1+\phi\right)+\Gamma\left[\bar{x}_{2}\left(\bar{x}_{1}+\bar{x}_{2}+\phi\left(\bar{x}_{12}+\bar{x}_{21}\right)\right)+\phi\left(\bar{s}_{d}-\bar{s}\right)\left(\bar{x}_{2}+\phi\bar{x}_{12}\right)\right]\right)}{\sigma_{0}\left(1+\phi\right)+\Gamma\left(\bar{x}_{1}+\bar{x}_{2}+\phi\left(\bar{x}_{12}+\bar{x}_{21}\right)\right)}$
(17b)

where $\Gamma=\frac{\beta_{0}\left(1+\phi\right)^{2}}{\phi^{2}}$. Although
these equations do not prove that the driven system synchronizes to the driver
system, they do indicate that the driven system rapidly approaches a manifold
in which the driven primary infectives are strongly affected by the infectives
in the driver system. We show in subsequent numerical simulations that the
driven variables do in fact approach those in the driver system.

##  4 Results

###  4.1 Comparison of Solutions

The original equations in the shifted, barred variables (Eqs. (21a)-(21j)) are
solved, and the values of the variables found on the right hand side of Eqs.
(16a)-(17b) are substituted into Eqs. (16a)-(17b) to find the value of
$\bar{x}_{1}$, $\bar{x}_{2}$, $\bar{x}_{d1}$, and $\bar{x}_{d2}$ given by the
center manifold equations. After computing in the shifted variables, we shift
back to the original variables.

Figure 2: Plot of $x_{1}$ computed using the complete system of equations
(blue) and the center manifold equation (red) for (a) $t=0$ to $t=100$, and
(b) $t=80$ to $t=100$. Parameter values are given by $\mu=0.02$,
$\sigma=50.0$, $\phi=3.0$, and $\beta=200.0$.

Figure 2(a) shows a comparison of $x_{1}$ found through computation of the
complete system and found through computation using the center manifold
equation. To show the agreement more clearly, Fig. 2(b) shows a portion of
Fig. 2(a). The prediction for the primary infectives is occasionally negative
since the center manifold equations involve $\bar{x}_{i}$ terms which
represent infectives shifted from their fixed point. As discussed in Shaw et
al. (2007), the predicted deviations may be large enough so that adding the
fixed point $x_{i,0}$ results in a negative value for the $x_{i}$ terms.
However, since the errors occur in regions where there is no epidemic
outbreak, they are not of biological importance. Furthermore, one can see that
the timing of the outbreaks predicted by the reduced system agrees well with
the actual outbreak time of the complete system.

Figure 3(a) and (c) shows the comparison of $x_{d1}$ and $x_{d2}$ respectively
found through computation of the complete system and found through computation
using the center manifold equation. The agreement is perfect and is seen more
clearly in Figs. 3(b) and (d) which show portion of Figs. 3(a) and (c)
respectively.

Figure 3: Plots of (a)-(b) $x_{d1}$ and (c)-(d) $x_{d2}$ computed using the
complete system of equations (blue) and the center manifold equation (red) for
(a) and (c) $t=0$ to $t=100$, and (b) and (d) $t=65$ to $t=100$. Parameter
values are given by $\mu=0.02$, $\sigma=50.0$, $\phi=3.0$, and $\beta=200.0$.

###  4.2 Transverse Stability

To show that the driven system synchronizes with the full system for the
components chosen, we can examine the transverse stability of the solution of
the subsystem. That is, we examine the difference between the full vector
field (the driver) and the subsystem. By computing the linear variational
equations of the difference, we can find the conditional Lyapunov exponents
for solutions near the subsystem.

To establish notation, we rewrite the components of the governing equations
for the two serotype multistrain disease model as

${\bf X}=[s,x_{1},x_{2}]^{T}\quad{\bf Y}=[r_{1},r_{2},x_{21},x_{12}]^{T},$
(18)

where $T$ denotes the transpose. Additionally, we can rewrite the components
of the governing equations for the two serotype multistrain disease subsystem
that are driven by the secondary infectious individuals as

${\bf Z}=[s_{d},x_{1d},x_{2d}]^{T}.$ (19)

To write the system of equations we consider the split differential equations

$\displaystyle{\bf{\bf\dot{X}}=F_{1}(X,Y)}$ (20a)
$\displaystyle{\bf{\bf\dot{Y}}=F_{2}(X,Y).}$ (20b) along with the subsystem
${\bf{\bf\dot{Z}}=F_{1}(Z,Y)}.$ (20c)

We would like ${\bf Z}(t)\rightarrow{\bf X}(t)$ in the asymptotic limit.
Letting ${\bm{\xi}}={\bf X}-{\bf Z}$, one obtains the linear variation about
${\bf X}(t)$ given by ${\dot{\bm{\xi}}}={\bf F}_{1,x}({\bf X},{\bf
Y}){\bm{\xi}}$. The Lyapunov exponents of the subsystem are the conditional
Lyapunov exponents (CLE). Negative CLE provide a sufficient condition for the
subsystem to converge, and are shown in Fig. 4.

Figure 4:  Time series of the conditional Lyapunov exponents. Parameter values
are given by $\mu=0.02$, $\sigma=50.0$, $\phi=3.0$, and $\beta=200.0$. The
initial condition for the driver system was a point on the chaotic attractor
obtained by computing a long time series, and the initial condition for the
driven system was a perturbation from the driver. Figure 5: Absolute error in
log variables for susceptibles and primary infectives versus time. Solid
green: $|\log s_{d}-\log s|$, solid blue: $|\log x_{1d}-\log x_{1}|$, dashed
red: $|\log x_{2d}-\log x_{2}|$. The initial condition for the driver system
was a point on the chaotic attractor obtained by computing a long time series.
The initial condition for the driven system was a random perturbation with
mean 10% away from the driver system. Parameter values are the same as in Fig.
4.

In Fig. 5, we show a sample time series of the error between the driven and
driver systems. Susceptibles and primary infectives in the driven system
rapidly converge to the values found in the driver system. Even when the
driven system is initially perturbed farther from the driver system,
convergence occurs within a similar time period.

##  5 Conclusions

In this article, we considered a two serotype multistrain disease model in
which susceptibles may develop a primary infection from either of the two
serotypes and upon recovery become immune to the strain that caused the
primary infection. However, the second serotype may cause a secondary
infection. We also consider a two serotype multistrain disease subsystem that
is driven by the secondary infectious individuals of the first model. By
performing center manifold analysis, we can reduce the dimension of these
driven and driver systems. In particular, we were able to analytically find an
approximate invariant manifold for the driven primary infectives expressed
almost entirely in terms of driver system variables. Numerical simulations
demonstrate the excellent agreement between solutions of the original, higher-
dimensional system and the lower-dimensional center manifold equations. We
further demonstrated the driven system is synchronized with the complete
system by examining the transverse stability of the subsystem solution.

In summary, we have developed a new method to determine unobserved epidemic
sub-populations using synchronization properties of the epidemic model. As an
application, it is known that the majority of dengue fever hospital cases are
secondary infections Nisalak et al. (2003). Thus one may have information
about secondary infections but not about primary infections. In driving a
model using the known secondary infection data, the primary infective
populations in the driven system synchronize to the correct values from the
driver system. It is therefore possible to deduce unobserved primary infection
levels in the population. It should be noted that this method probably works
for other multistrain models, such as the dengue model in Bianco et al. (2009)
in which we include temporary cross immunity immediately following a primary
infection.

Future work includes the application of our method to real disease data.
Before this can be performed, one would first have to perform the analysis in
the presence of noise (i.e. the driver signal comes from a noisy system and/or
has additional errors in it). Additionally, the theory will need to be applied
to the situation when data is sampled discretely. Preliminary testing using
discretely sampled data for the secondary infectives has shown that the
technique can indeed recover primary infectives.

###### Acknowledgements.

EF is supported by Award Number CMMI-1233397 from the National Science
Foundation. LBS is supported by Award Number R01GM090204 from the National
Institute of General Medical Sciences. The content is solely the
responsibility of the authors and does not necessarily represent the official
views of the National Institute Of General Medical Sciences or the National
Institutes of Health. IBS is supported by the NRL Base Research Program
contract number N0001414WX00023 and by the Office of Naval Research contract
number N0001414WX20610.

## References

  * Anderson and May (1991) Anderson, R. M. and R. M. May, 1991.  Infectious Diseases of Humans.  Oxford University Press.
  * Bianco et al. (2009) Bianco, S., L. B. Shaw and I. B. Schwartz, 2009.  Epidemics with multistrain interactions: The interplay between cross immunity and antibody-dependent enhancement.  Chaos 19, 043123.
  * Billings et al. (2007) Billings, L. et al., 2007.  Instabilities in multiserotype disease models with antibody-dependent enhancement.  J. Theor. Biol. 246, 18.
  * Bjornstad et al. (2002) Bjornstad, O. N., B. F. Finkenstadt and B. T. Grenfell, 2002.  Dynamics of measles epidemics: Estimating scaling of transmission rates using a time series sir model.  Ecological Monographs 72(2), 169–184.
  * Blarer and Doebeli (1999) Blarer, A. and M. Doebeli, 1999.  Resonance effects and outbreaks in ecological time series.  Ecol. Lett. 2, 167–177.
  * Boccaletti et al. (2002) Boccaletti, S. et al., 2002.  The synchronization of chaotic systems.  Physics Reports-Review Section Of Physics Letters 366, 1–101.
  * Carr (1981) Carr, J., 1981.  Applications of Centre Manifold Theory.  Springer-Verlag.
  * Chen and Lü (2002) Chen, S. and J. Lü, 2002.  Parameters identification and synchronization of chaotic systems based upon adaptive control.  Physics Letters A 299, 353.
  * Chicone and Latushkin (1997) Chicone, C. and Y. Latushkin, 1997.  Center manifolds for infinite dimensional nonautonomous differential equations.  J. Differ. Equations 141, 356–399.
  * Cummings et al. (2005) Cummings, D. A. T. et al., 2005.  Dynamic effects of antibody-dependent enhancement on the fitness of viruses.  Proc. Natl. Acad. Sci. U.S.A. 102(42), 15259–15264.
  * Dedieu and Ogorzalek (1997) Dedieu, H. and M. J. Ogorzalek, 1997.  Identifiability and identification of chaotic systems based on adaptive synchronization.  IEEE transactions on circuits and systems—I: Fundamental theory and applications 44, 948.
  * Duan et al. (2003) Duan, J., K. Lu and B. Schmalfuss, 2003.  Inavriant manifolds for stochastic partial differential equations.  Ann. Probab. 31(4), 2109–2135.
  * Ferguson et al. (1999) Ferguson, N. M., C. A. Donnelly and R. M. Anderson, 1999.  Transmission dynamics and epidemiology of dengue: insights from age-stratified sero-prevalence surveys.  Philos. Trans. R. Soc. Lond. Ser. B, Biol. Sci. 354, 757–768.
  * Forgoston et al. (2009) Forgoston, E., L. Billings and I. B. Schwartz, 2009.  Accurate noise projection for reduced stochastic epidemic models.  Chaos 19, 043110.
  * Forgoston and Schwartz (2013) Forgoston, E. and I. B. Schwartz, 2013.  Predicting unobserved exposures from seasonal epidemic data.  Bull Math Biol 75, 1450.
  * Gibson et al. (2004) Gibson, G. J., A. Kleczkowski and C. A. Gilligan, 2004.  Bayesian analysis of botanical epidemics using stochastic compartmental models.  PNAS 101, 12120.
  * Huang and Lin (2013) Huang, L. and L. Lin, 2013.  Parameter identification and synchronization of uncertain chaotic systems based on sliding mode observer.  Mathematical Problems in Engineering 2013, 859304.
  * Lekone and Finkenstädt (2006) Lekone, P. E. and B. F. Finkenstädt, 2006.  Statistical inference in a stochastic epidemic seir model with control intervention: Ebola as a case study.  Biometrics 62, 1170.
  * Nagao and Koelle (2008) Nagao, Y. and K. Koelle, 2008.  Decreases in dengue transmission may act to increase the incidence of dengue hemorrhagic fever.  Proc. Nat. Acad. Sci. USA 105, 2238–2243.
  * Nisalak et al. (2003) Nisalak, A. et al., 2003.  Serotype-specific dengue virus circulation and dengue disease in Bangkok, Thailand from 1973 to 1999.  Am. J. Trop. Med. Hyg. 68, 191–202.
  * Parlitz et al. (1996) Parlitz, U., L. Junge and L. Kocarev, 1996.  Synchronization-based parameter estimation from time series.  Physical Review E 54, 6253.
  * Rigau-Perez et al. (1998) Rigau-Perez, J. et al., 1998.  Dengue and dengue hemmorrhagic fever.  Lancet 352, 971–977.
  * Roberts (2008) Roberts, A. J., 2008.  Normal form transforms separate slow and fast modes in stochastic dynamical systems.  Physica A 387(1), 12–38.
  * Schaffer et al. (1993) Schaffer, W. M. et al., 1993.  Transient periodicity and episodic predictability in biological dynamics.  IMA J. Math. Appl. Med. 10, 227–247.
  * Schwartz and Smith (1983) Schwartz, I. and H. Smith, 1983.  Infinite subharmonic bifurcations in an SEIR epidemic model.  J. Math. Biol. 18, 233–253.
  * Schwartz et al. (2005) Schwartz, I. B. et al., 2005.  Chaotic desynchronization of multi-strain diseases.  Phys. Rev. E 72, 066201.
  * Shaw et al. (2007) Shaw, L. B., L. Billings and I. B. Schwartz, 2007.  Using dimension reduction to improve outbreak predictability of multistrain diseases.  J. Math. Bio. 55, 1–19.

##  Appendix A Shifted, Rescaled, and Augmented System of Equations

The governing equations for the two serotype multistrain disease model and the
subsystem driven by secondary infectious individuals are given by Eqs.
(4a)-(5c). We define a new set of variables, $\bar{s}$, $\bar{x}_{i}$,
$\bar{r}_{i}$, $\bar{x}_{ij}$, $\bar{s}_{d}$, and $\bar{x}_{id}$ for all $i,j$
as $\bar{s}(t)=s(t)-s_{0}$, $\bar{x}_{i}(t)=x_{i}(t)-x_{i,0}$,
$\bar{r}_{i}(t)=r(t)-r_{i,0}$, $\bar{x}_{ij}(t)=x_{ij}(t)-x_{ij,0}$,
$\bar{s}_{d}(t)=s_{d}(t)-s_{d,0}$, $\bar{x}_{id}(t)=x_{id}(t)-x_{id,0}$, and
these new variables are substituted into Eqs. (4a)-(5c).

Then, treating $\mu$ as a small parameter, we rescale time by letting
$t=\mu\tau$. We may then introduce the following rescaled parameters:
$\beta=\beta_{0}/\mu$ and $\sigma=\sigma_{0}/\mu$, where $\beta_{0}$ and
$\sigma_{0}$ are $\mathcal{O}(1)$. The inclusion of the parameter $\mu$ as a
new state variable means that the terms in our rescaled system which contain
$\mu$ are now nonlinear terms. Furthermore, the system is augmented with the
auxiliary equation $\frac{d\mu}{d\tau}=0$. The addition of this auxiliary
equation contributes an extra simple zero eigenvalue to the system and adds
one new center direction that has trivial dynamics. The shifted and rescaled,
augmented system of equations is given as

$\displaystyle\frac{d\bar{s}}{d\tau}=\mu^{2}-\beta_{0}\left(\bar{s}+\frac{\sigma_{0}}{\beta_{0}(1+\phi)}\right)\left(\bar{x}_{1}+\bar{x}_{2}+\phi\left(\bar{x}_{21}+\bar{x}_{12}\right)+\frac{\mu^{2}(1+\phi)}{\sigma_{0}}\right),$
(21a)
$\displaystyle\frac{d\bar{x}_{1}}{d\tau}=\beta_{0}\left(\bar{s}+\frac{\sigma_{0}}{\beta_{0}(1+\phi)}\right)\left(\bar{x}_{1}+\phi\bar{x}_{21}+\frac{\mu^{2}(1+\phi)}{2\sigma_{0}}\right)-\sigma_{0}\left(\bar{x}_{1}+\frac{\mu^{2}}{2\sigma_{0}}\right),$
(21b)
$\displaystyle\frac{d\bar{x}_{2}}{d\tau}=\beta_{0}\left(\bar{s}+\frac{\sigma_{0}}{\beta_{0}(1+\phi)}\right)\left(\bar{x}_{2}+\phi\bar{x}_{12}+\frac{\mu^{2}(1+\phi)}{2\sigma_{0}}\right)-\sigma_{0}\left(\bar{x}_{2}+\frac{\mu^{2}}{2\sigma_{0}}\right),$
(21c)
$\displaystyle\frac{d\bar{r}_{1}}{d\tau}=\sigma_{0}\left(\bar{x}_{1}+\frac{\mu^{2}}{2\sigma_{0}}\right)-\beta_{0}\left(\bar{r}_{1}+\frac{\sigma_{0}}{\beta_{0}(1+\phi)}\right)\left(\bar{x}_{2}+\phi\bar{x}_{12}+\frac{\mu^{2}(1+\phi)}{2\sigma_{0}}\right),$
(21d)
$\displaystyle\frac{d\bar{r}_{2}}{d\tau}=\sigma_{0}\left(\bar{x}_{2}+\frac{\mu^{2}}{2\sigma_{0}}\right)-\beta_{0}\left(\bar{r}_{2}+\frac{\sigma_{0}}{\beta_{0}(1+\phi)}\right)\left(\bar{x}_{1}+\phi\bar{x}_{21}+\frac{\mu^{2}(1+\phi)}{2\sigma_{0}}\right),$
(21e)
$\displaystyle\frac{d\bar{x}_{21}}{d\tau}=\beta_{0}\left(\bar{r}_{2}+\frac{\sigma_{0}}{\beta_{0}(1+\phi)}\right)\left(\bar{x}_{1}+\phi\bar{x}_{21}+\frac{\mu^{2}(1+\phi)}{2\sigma_{0}}\right)-\sigma_{0}\left(\bar{x}_{21}+\frac{\mu^{2}}{2\sigma_{0}}\right),$
(21f)
$\displaystyle\frac{d\bar{x}_{12}}{d\tau}=\beta_{0}\left(\bar{r}_{1}+\frac{\sigma_{0}}{\beta_{0}(1+\phi)}\right)\left(\bar{x}_{2}+\phi\bar{x}_{12}+\frac{\mu^{2}(1+\phi)}{2\sigma_{0}}\right)-\sigma_{0}\left(\bar{x}_{12}+\frac{\mu^{2}}{2\sigma_{0}}\right),$
(21g)
$\displaystyle\frac{d\bar{s}_{d}}{d\tau}=\mu^{2}-\beta_{0}\left(\bar{s}_{d}+\frac{\sigma_{0}}{\beta_{0}(1+\phi)}\right)\left(\bar{x}_{1d}+\bar{x}_{2d}+\phi\left(\bar{x}_{21}+\bar{x}_{12}\right)+\frac{\mu^{2}(1+\phi)}{\sigma_{0}}\right),$
(21h)
$\displaystyle\frac{d\bar{x}_{1d}}{d\tau}=\beta_{0}\left(\bar{s}_{d}+\frac{\sigma_{0}}{\beta_{0}(1+\phi)}\right)\left(\bar{x}_{1d}+\phi\bar{x}_{21}+\frac{\mu^{2}(1+\phi)}{2\sigma_{0}}\right)-\sigma_{0}\left(\bar{x}_{1d}+\frac{\mu^{2}}{2\sigma_{0}}\right),$
(21i)
$\displaystyle\frac{d\bar{x}_{2d}}{d\tau}=\beta_{0}\left(\bar{s}_{d}+\frac{\sigma_{0}}{\beta_{0}(1+\phi)}\right)\left(\bar{x}_{2d}+\phi\bar{x}_{12}+\frac{\mu^{2}(1+\phi)}{2\sigma_{0}}\right)-\sigma_{0}\left(\bar{x}_{2d}+\frac{\mu^{2}}{2\sigma_{0}}\right),$
(21j) $\displaystyle\frac{d\mu}{d\tau}=0,$ (21k)

where the endemic fixed point is now located at the origin.

##  Appendix B Definition of New Variables

Using the fact that
$(\bar{s},\bar{x}_{1},\bar{x}_{2},\bar{r}_{1},\bar{r}_{2},\bar{x}_{21},\bar{x}_{12},\bar{s}_{d},\bar{x}_{1d},\bar{x}_{2d})^{T}={\bf
P}\cdot{\bf W}^{T}$, where ${\bf P}$ is given by Eq. (7) and ${\bf
W}=(W_{1},W_{2},W_{3},W_{4},W_{5},W_{6},W_{7},W_{8},W_{9},W_{10})$, then the
transformation matrix leads to the following definition of new variables,
$W_{i}$, $i=1\ldots 10$:

$\begin{split}&W_{1}=\frac{\bar{x}_{21}-\bar{x}_{1}}{1+\phi},\quad
W_{2}=\frac{\bar{x}_{12}-\bar{x}_{2}}{1+\phi},\quad
W_{3}=\frac{\bar{x}_{1d}+\bar{x}_{2d}-\bar{x}_{1}-\bar{x}_{2}}{\phi},\quad
W_{4}=\bar{x}_{2d}-\bar{x}_{2},\\\ &W_{5}=\bar{s},\quad
W_{6}=\frac{\bar{x}_{1}+\phi\bar{x}_{21}}{1+\phi},\quad
W_{7}=\frac{\bar{x}_{2}+\phi\bar{x}_{12}}{1+\phi},\quad
W_{8}=\frac{\phi\bar{r}_{1}+\phi\bar{x}_{1}+\bar{r}_{1}-\phi\bar{x}_{21}}{1+\phi},\\\
&W_{9}=\frac{\phi\bar{r}_{2}+\phi\bar{x}_{2}+\bar{r}_{2}-\phi\bar{x}_{12}}{1+\phi},\quad
W_{10}=\frac{\phi\bar{s}_{d}-\bar{x}_{1d}-\bar{x}_{2d}+\bar{x}_{1}+\bar{x}_{2}}{\phi}.\end{split}$
(22)

##  Appendix C Transformed Evolution Equations

The application of the transformation matrix ${\bf P}$ given by Eq. (7) to
Eqs. (21a)-(21j) leads to the following set of transformed evolution
equations:

$\displaystyle\frac{dW_{1}}{d\tau}=\beta_{0}\left(W_{6}+\frac{\mu^{2}}{2\sigma_{0}}\right)\left(W_{9}+\phi
W_{2}-W_{5}\right)-\sigma_{0}W_{1},$ (23a)
$\displaystyle\frac{dW_{2}}{d\tau}=\beta_{0}\left(W_{7}+\frac{\mu^{2}}{2\sigma_{0}}\right)\left(W_{8}+\phi
W_{1}-W_{5}\right)-\sigma_{0}W_{2},$ (23b)
$\displaystyle\frac{dW_{3}}{d\tau}=W_{3}\left(\beta_{0}\left(W_{3}+W_{10}\right)-\frac{\sigma_{0}\phi}{1+\phi}\right)$
$\displaystyle\hskip
28.45274pt+\frac{\beta_{0}\left(1+\phi\right)}{\phi}\left(W_{6}+W_{7}+\frac{\mu^{2}}{\sigma_{0}}\right)\left(W_{3}+W_{10}-W_{5}\right),$
(23c)
$\displaystyle\frac{dW_{4}}{d\tau}=\left(\beta_{0}\left(W_{3}+W_{10}\right)-\frac{\sigma_{0}\phi}{1+\phi}\right)W_{4}$
$\displaystyle\hskip
28.45274pt+\beta_{0}\left(W_{3}+W_{10}-W_{5}\right)\left(\left(1+\phi\right)W_{7}+\frac{\mu^{2}\left(1+\phi\right)}{2\sigma_{0}}\right),$
(23d)
$\displaystyle\frac{dW_{5}}{d\tau}=\mu^{2}-\beta_{0}\left(W_{5}+\frac{\sigma_{0}}{\beta_{0}(1+\phi)}\right)\left((1+\phi)\left(W_{6}+W_{7}+\frac{\mu^{2}}{\sigma_{0}}\right)\right),$
(23e)
$\displaystyle\frac{dW_{6}}{d\tau}=\beta_{0}\left(W_{6}+\frac{\mu^{2}}{2\sigma_{0}}\right)\left(\phi^{2}W_{2}+\phi
W_{9}+W_{5}\right),$ (23f)
$\displaystyle\frac{dW_{7}}{d\tau}=\beta_{0}\left(W_{7}+\frac{\mu^{2}}{2\sigma_{0}}\right)\left(\phi^{2}W_{1}+\phi
W_{8}+W_{5}\right),$ (23g)
$\displaystyle\frac{dW_{8}}{d\tau}=\left(W_{6}+\frac{\mu^{2}}{2\sigma_{0}}\right)\left(\sigma_{0}+\beta_{0}\left(\phi
W_{5}-\phi W_{9}-\phi^{2}W_{2}\right)\right)$ $\displaystyle\hskip
28.45274pt-\left(1+\phi\right)\beta_{0}\left(W_{7}+\frac{\mu^{2}}{2\sigma_{0}}\right)\left(\phi
W_{1}+W_{8}+\frac{\sigma_{0}}{\beta_{0}\left(1+\phi\right)}\right),$ (23h)
$\displaystyle\frac{dW_{9}}{d\tau}=\left(W_{7}+\frac{\mu^{2}}{2\sigma_{0}}\right)\left(\sigma_{0}+\beta_{0}\left(\phi
W_{5}-\phi W_{8}-\phi^{2}W_{1}\right)\right)$ $\displaystyle\hskip
28.45274pt-\left(1+\phi\right)\beta_{0}\left(W_{6}+\frac{\mu^{2}}{2\sigma_{0}}\right)\left(\phi
W_{2}+W_{9}+\frac{\sigma_{0}}{\beta_{0}\left(1+\phi\right)}\right),$ (23i)
$\displaystyle\frac{dW_{10}}{d\tau}=\mu^{2}-\left(1+\phi\right)\beta_{0}W_{3}\left(W_{3}+W_{10}\right)$
$\displaystyle\hskip
28.45274pt+\frac{\left(1+\phi\right)\beta_{0}}{\phi}\left(W_{6}+W_{7}+\frac{\mu^{2}}{\sigma_{0}}\right)\left(W_{5}-\left(1+\phi\right)\left(W_{3}+W_{10}\right)-\frac{\phi\sigma_{0}}{\beta_{0}\left(1+\phi\right)}\right),$
(23j) $\displaystyle\frac{d\mu}{d\tau}=0.$ (23k)

##  Appendix D Center Manifold Condition

Substitution of the center manifold functions $W_{i}=h_{i}$ given by Eq. (12)
into the transformed evolution equations given in Appendix C leads to the
following center manifold condition:

$\displaystyle\frac{\partial h_{1}}{\partial
W_{5}}\frac{dW_{5}}{d\tau}+\frac{\partial h_{1}}{\partial
W_{6}}\frac{dW_{6}}{d\tau}+\frac{\partial h_{1}}{\partial
W_{7}}\frac{dW_{7}}{d\tau}+\frac{\partial h_{1}}{\partial
W_{8}}\frac{dW_{8}}{d\tau}+\frac{\partial h_{1}}{\partial
W_{9}}\frac{dW_{9}}{d\tau}+\frac{\partial h_{1}}{\partial
W_{10}}\frac{dW_{10}}{d\tau}$ $\displaystyle\hskip
28.45274pt=\beta_{0}\left(W_{6}+\frac{\mu^{2}}{2\sigma_{0}}\right)\left(W_{9}+\phi
h_{2}-W_{5}\right)-\sigma_{0}h_{1},$ (24a) $\displaystyle\frac{\partial
h_{2}}{\partial W_{5}}\frac{dW_{5}}{d\tau}+\frac{\partial h_{2}}{\partial
W_{6}}\frac{dW_{6}}{d\tau}+\frac{\partial h_{2}}{\partial
W_{7}}\frac{dW_{7}}{d\tau}+\frac{\partial h_{2}}{\partial
W_{8}}\frac{dW_{8}}{d\tau}+\frac{\partial h_{2}}{\partial
W_{9}}\frac{dW_{9}}{d\tau}+\frac{\partial h_{2}}{\partial
W_{10}}\frac{dW_{10}}{d\tau}$ $\displaystyle\hskip
28.45274pt=\beta_{0}\left(W_{7}+\frac{\mu^{2}}{2\sigma_{0}}\right)\left(W_{8}+\phi
h_{1}-W_{5}\right)-\sigma_{0}h_{2},$ (24b) $\displaystyle\frac{\partial
h_{3}}{\partial W_{5}}\frac{dW_{5}}{d\tau}+\frac{\partial h_{3}}{\partial
W_{6}}\frac{dW_{6}}{d\tau}+\frac{\partial h_{3}}{\partial
W_{7}}\frac{dW_{7}}{d\tau}+\frac{\partial h_{3}}{\partial
W_{8}}\frac{dW_{8}}{d\tau}+\frac{\partial h_{3}}{\partial
W_{9}}\frac{dW_{9}}{d\tau}+\frac{\partial h_{3}}{\partial
W_{10}}\frac{dW_{10}}{d\tau}$ $\displaystyle\hskip
28.45274pt=h_{3}\left(\beta_{0}\left(h_{3}+W_{10}\right)-\frac{\sigma_{0}\phi}{1+\phi}\right)$
$\displaystyle\hskip
28.45274pt+\frac{\beta_{0}\left(1+\phi\right)}{\phi}\left(W_{6}+W_{7}+\frac{\mu^{2}}{\sigma_{0}}\right)\left(h_{3}+W_{10}-W_{5}\right),$
(24c) $\displaystyle\frac{\partial h_{4}}{\partial
W_{5}}\frac{dW_{5}}{d\tau}+\frac{\partial h_{4}}{\partial
W_{6}}\frac{dW_{6}}{d\tau}+\frac{\partial h_{4}}{\partial
W_{7}}\frac{dW_{7}}{d\tau}+\frac{\partial h_{4}}{\partial
W_{8}}\frac{dW_{8}}{d\tau}+\frac{\partial h_{4}}{\partial
W_{9}}\frac{dW_{9}}{d\tau}+\frac{\partial h_{4}}{\partial
W_{10}}\frac{dW_{10}}{d\tau}$ $\displaystyle\hskip
28.45274pt=\left(\beta_{0}\left(h_{3}+W_{10}\right)-\frac{\sigma_{0}\phi}{1+\phi}\right)h_{4}$
$\displaystyle\hskip
28.45274pt+\beta_{0}\left(h_{3}+W_{10}-W_{5}\right)\left(\left(1+\phi\right)W_{7}+\frac{\mu^{2}\left(1+\phi\right)}{2\sigma_{0}}\right).$
(24d)