Identifiability and Observability Analysis for Epidemiological Models: Insights on the SIRS Model

Alicja B. Kubik akubik@ucm.es MOMAT, Instituto de Matemática Interdisciplinar (IMI), Univ. Complutense de Madrid, 28040 Madrid, Spain Benjamin Ivorra ivorra@mat.ucm.es MOMAT, Instituto de Matemática Interdisciplinar (IMI), Univ. Complutense de Madrid, 28040 Madrid, Spain Alain Rapaport alain.rapaport@inrae.fr MISTEA, Univ. Montpellier, INRAE, Institut Agro, 34060 Montpellier, France Ángel M. Ramos angel@mat.ucm.es MOMAT, Instituto de Matemática Interdisciplinar (IMI), Univ. Complutense de Madrid, 28040 Madrid, Spain

Abstract

The problems of observability and identifiability have been of great interest as previous steps to estimating parameters and initial conditions of dynamical systems to which some known data (observations) are associated. While most works focus on linear and polynomial/rational systems of ODEs, general nonlinear systems have received far less attention and, to the best of our knowledge, no general constructive methodology has been proposed to assess and guarantee parameter and state recoverability in this context. We consider a class of systems of parameterized nonlinear ODEs and some observations, and study if a system of this class is observable, identifiable or jointly observable-identifiable; our goal is to identify its parameters and/or reconstruct the initial condition from the data. To achieve this, we introduce a family of efficient and fully constructive procedures that allow recoverability of the unknowns with low computational cost and address the aforementioned gap. Each procedure is tailored to different observational scenarios and based on the resolution of linear systems. As a case study, we apply these procedures to several epidemic models, with a detailed focus on the SIRS model, demonstrating its joined observability-identifiability when only a portion of the infected individuals is measured, a scenario that has not been studied before. In contrast, for the same observations, the SIR model is observable and identifiable, but not jointly observable-identifiable. This distinction allows us to introduce a novel approach to discriminating between different epidemiological models (SIR vs. SIRS) from short-time data. For these two models, we illustrate the theoretical results through some numerical experiments, validating the approach and highlighting its practical applicability to real-world scenarios.

Key words: Nonlinear epidemiological models; Parameter identifiability; State observability; Joined observability-identifiability; SIRS model analysis; Epidemiological model discrimination; Epidemic data reconstruction.

1 Introduction

Mathematical modeling has been widely used for studying epidemics for a long time. Some of the most popular models in epidemics generated by infectious diseases are compartmental models, whose theory was set up in the early 1900s (see [40]), and have been regularly used since then to estimate the dynamics of different diseases (for example, influenza [7, 11], tuberculosis [20, 36], Ebola [17, 61], or vector-transmitted diseases such as malaria [62], dengue fever or the Zika virus [8, 21]), with a significant surge in applications in recent years due to the COVID-19 pandemic (see, for example, [2, 6, 33, 54, 56]). These models are usually deterministic systems of nonlinear ODEs, which include the class of systems that we are going to address from now on.

Given a disease and some information on the different biological processes that are involved, we need to set up a model that captures the most important features of its dynamics, being careful not to complexify it up to a point that may result intractable from a mathematical perspective. To these dynamical models there are associated data, available in the literature (e.g., parameter values) or collected over time (e.g., how many people are hospitalized or vaccinated at certain dates). This information is partial most of the time, in the sense that they are related to some of the parameters and state variables of the model, but rarely to all of them. For instance, the sizes of susceptible or asymptomatic compartments are often hard or even impossible to measure (they are hidden variables), and the transmission rate of a new virus is often not known in advance.

Once we have settled a model and know some information about the parameters and compartments, a question arises: Can we uniquely reconstruct all the unknowns of the model from partial observations? The ability to accurately reconstruct key aspects of disease dynamics from observed data is fundamental to understand epidemic trajectories and design effective control strategies. Thus, determining whether the parameters can be uniquely identified or the states can be precisely reconstructed from available data is essential. We thus face an inverse problem. Such identification and state reconstruction problems are often met in automatic control for linear and nonlinear dynamics. Although many successful applications have been recorded in aeronautics, automotive, electronics, or chemistry, among other domains, very few has been investigated for epidemiological models. Moreover, despite substantial research addressing these questions, significant challenges remain, particularly in nonlinear models, where exact reconstruction of parameters and states is often elusive and requires sophisticated analytical approaches.

In this context, the determination of a feasible set of parameters and an initial condition that align with observed data falls into two key problems: The parameter identification problem and the state observation problem in automatic control literature. These have been addressed in many different ways (see, for example, [10, 12, 29, 14, 60, 51, 64]). However, before attempting to determine these unknowns, it is crucial to assess the extent of recoverable information from the available epidemic data. For example, can we determine the disease contact rate if we know the data of hospitalized people? What about the loss of immunity rate? How many people were infected when the data started to be reported? Of course, these are not issues related only to epidemiological models. In general, given a phenomenon we want to model (e.g., physical, biological, or mechanical), we will have some parameters, a state vector (in the epidemiological case, these states are the size of the population in each compartment, or a portion of this size), and some observed data of this phenomenon (also called measurements, outputs of the model, or signals [38]). The theories of identifiability and observability provide the foundation for addressing these questions. Let us briefly introduce these concepts.

Identifiability theory is the one in charge of deciding, given an initial condition and some observed data, if we can determine univocally all the unknown parameters that govern our system. If it is not possible, we can possibly recover them partially (some of them or some combination of them).

Observability theory is the one in charge of deciding, given the parameters and some observed data, if we can recover the initial condition; in other words, if there exists a unique initial condition such that the solution starting from this initial condition matches our observations, and hence we can distinguish different states from partial measurements of the state vector.

If one is interested in both identifying and observing the system, i.e., if it is jointly observable-identifiable (see [13]) it is common to extend the system by considering the parameters as part of the state vector (with a null dynamics), and then considering the observation problem of this extended system. The most commonly used technique to decide whether a system is observable or not is the Hermann-Krener condition, or Observability Rank Condition, that consists in studying the separability of a set composed by the observations and its Lie derivatives with respect to the vector field of the system (see [31]); this computation may be sometimes simplified by exploiting symmetries or groups of invariance of the system (in particular, for mechanical and/or robotics systems, see e.g. [48]). If we conclude our system is observable, then we can make use of several techniques that may help us estimate practically the unknowns; in particular, to perform this task, we can try to construct a state observer, as for example the Luenberger observer [47], high-gains observers [24], or the Kalman filter [37]. These techniques aim to estimate the state vector at a certain speed with a certain accuracy, minimizing the estimation error at a given time, such that it is usually asymptotically null, but do not guarantee exact reconstruction. Sometimes, however, one can treat the system algebraically and try to reconstruct exactly the unknowns in terms of, for example, the derivatives of the data whenever they exist and are known (or can be computed) perfectly (see [13, Chapter 3]).

The theories of identifiability and observability extend beyond the epidemiological context and are relevant to any phenomena modeled through ODE systems with measurable outputs or functions of the states (see [3, 4, 25, 32, 34, 65]). Indeed, there is few literature about studying the joined observability-identifiability of epidemiological models (see [13], which is a survey that can be used as a textbook about this topic, and [19], [30], and [49]).

Several works have emphasized the importance of initial conditions for identifying the parameters of a general system, and not every initial condition is suitable, i.e., some initial conditions can produce the same output when considering different parameters (see, e.g., [15], [59]). While there are fields where one can perform different experiments with chosen initial conditions in order to identify the parameters (see, e.g., [15, 35]), in disciplines such as Epidemiology this is not possible, and hence we need to study which initial conditions are not suitable, and if we can avoid them. In this paper, we build on the theory developed in [44], and extend the results by considering the existence of non-suitable initial conditions which are dependent on the parameters.

On the other hand, most of the research carried out in observability and identifiability of nonlinear dynamics is centered on polynomial and rational systems of equations (see, e.g., [16, 26, 28, 27, 59]), and the study of nonlinear, non-rational systems is scarce (see the survey [1]). In rational systems, the goal is typically finding the input-output (IO) equations in order to determine if the parameters are identifiable. These are structural equations of the system under study, which relate the parameters, the outputs and their derivatives, and their inputs (if any). However, we consider that, in these works, there lacks a direct, unified argument to conclude that these parameters can be effectively recovered from these equations. In [44], general nonlinearities are considered, and it advances the state of the art by minimizing reliance on high-order derivatives, improving robustness with respect to noisy and sparse datasets, and offering novel techniques to address joined observability-identifiability in nonlinear systems of ODEs. Moreover, this lacking argument is also covered by considering an assumption about linear independence of functions of the observations and their derivatives. This assumption of linear independence was already tackled in [15] to prove identifiability alone, however, it is framed in the case of rational systems and it assumes knowledge of initial conditions and having data at the initial time; these last assumptions are not feasible in some fields such as Epidemiology. It is as well mentioned in [53] for linear ODE systems, which do not cover epidemiological models either.

This paper is structured as follows. In Section 2, we motivate through the SIRS model the application of this framework to epidemiological models. Then, in Section 3, we revisit and extend the theoretical part presented in [44], and establish different constructive algorithms to recover the parameters and/or initial condition. We also extend the theory to incorporate higher-order derivatives, and highlight that this framework can be extended to scenarios with piecewise constant parameters, provided that the time instants of parameter changes are known or the time intervals with constant values are large enough.

In Section 4, we validate this approach by applying it to the SIRS model, proving joined observability-identifiability under partial observations of the infected compartment, an example that we have not seen considered in the literature. We will observe that in this case it is crucial considering initial conditions in parameter-dependent sets. Furthermore, although this is a rational model which can be studied through classical techniques, using our methodology allows for early-stage model discrimination between SIR and SIRS structures, enabling the detection of immunity loss directly from observational data. This distinction addresses a key gap in the literature, to the best of our knowledge, where theoretical studies often overlook model selection based on limited data. After this, we present in Section 5 other epidemiological examples to highlight the interest of the proposed approach. Finally, we perform in Section 6 different numerical experiments using the SIRS and the SIR models in two different scenarios that illustrate our approach.

2 Motivating framework: Insights from the SIRS model

Compartmental models, such as the SIRS model, strongly depend on real data to determine their parameters and the current state of the disease. These data are typically scarce and noisy, making it crucial to develop methods to estimate the unknown parameters and/or initial condition from partial observations. For this, we are going to follow two different approaches, depending on the type of available data: identifiability and observability, both formally defined in Section 3.

Let us consider the following classical SIRS model together with an output, given by an unknown fraction $k\in(0,1]$ of the infected individuals at time $t\geq 0$ , represented by $y$ :

\left\{\begin{array}[]{ccl}\dot{S}&=&-\beta SI+\mu R,\\ \dot{I}&=&\beta SI-\gamma I,\\ \dot{R}&=&\gamma I-\mu R,\\ y&=&kI,\end{array}\quad\forall\,t\geq 0.\right.

(1)

This model is particularly relevant for understanding the dynamics of diseases such as influenza of some coronaviruses. It is also of interest in different realistic situations, such as performing random tests among the population, or considering that a fraction $1-k$ of the infectious people are asymptomatic and hence we do not measure those cases, where the parameter $k$ is typically unknown. For instance, during the COVID-19 pandemic, one of the main challenges when studying the prevalence of the disease was estimating the asymptomatic cases, whereas governments also conducted randomized testing to estimate this prevalence (see, for example, [55]).

Given this model, we are interested in knowing whether it is identifiable, observable or jointly observable-identifiable. To address this, in Section 3, we revisit and extend the general framework of [44]. The presented SIRS model with output $kI$ fits within the class of systems presented in Section 3 and serves as a motivating example. Indeed, in Section 4 we make use of this approach to prove its joined observability-identifiability when $\mu>0$ . Moreover, we compare these results to those obtained for the (a priori simpler) SIR model, which can be considered as a limiting case of the SIRS model with $\mu=0$ , considering the same observational data (i.e., a fraction of infected individuals). Our analysis reveals that the SIR model is both observable and identifiable, but lacks joined observability-identifiability, highlighting a key distinction between these two models that allows us to distinguish them in an early stage. We will illustrate this difference through the numerical tests presented in Section 6.

3 A general framework

In this section, we extend the general framework established in [44] for analyzing observability, identifiability and joined observability-identifiability of dynamical systems, particularly focusing on systems modeled by autonomous differential equations with unknown parameters. We provide conditions under which the parameters and/or initial condition of such systems can be uniquely determined from observations of the system’s output. To this end, we first introduce the mathematical formulation of the model and provide formal definitions of the key concepts: Identifiability, observability, and joined observability-identifiability. Unlike in [44], when discussing identifiability, we will consider the possibility of parameter-dependent initial conditions that impede identifiability; this is a typical situation in epidemiological models. Additionally, we revisit some concepts on Lie derivatives that are instrumental in the theoretical analysis. Building on this foundation, we extend the key results from [44], taking into account these parameter-dependent initial conditions, that enable to establish different constructive methodologies to recover the initial condition and/or the parameters, whenever possible. Moreover, we present a third approach using higher-order derivatives to recover the unknowns in cases where the system is proven to be jointly observable-identifiable.

3.1 Mathematical formulation and definitions

Consider a phenomenon that can be modeled by a system of first order autonomous differential equations, which depends on some unknown parameter vector $\theta$ . Together with an initial condition $\xi$ , we assume the model provides some output $y_{(\xi,\theta)}$ , described by a suitable function $h$ . The mathematical formulation of the model and its output is given by

\left\{\begin{array}[]{l}\dfrac{\mathrm{d}x}{\mathrm{d}t}(t;\xi,\theta)=f(x(t;\xi,\theta),\theta),\quad x(0;\xi,\theta)=\xi,\\[10.00002pt] y_{(\xi,\theta)}(t)=h(x(t;\xi,\theta),\theta),\end{array}\right.

(2)

where $f(\cdot,\cdot):\Omega\times\Theta\rightarrow\mathbb{R}^{n}$ is a known function of $(x,\theta)$ which is locally Lipschitz-continuous w.r.t. $x\in\Omega\subset\mathbb{R}^{n}$ (to guarantee uniqueness of solutions) and continuous w.r.t. $\theta\in\Theta\subset\mathbb{R}^{b}$ ; $\theta$ are the constant parameters of the system; $\Omega\subset\mathbb{R}^{n}$ is a positively invariant set with respect to the system of ODEs of System (2), for any $\theta\in\Theta$ (i.e., solutions starting in $\Omega$ remain in $\Omega$ during all its definition time); $x(\cdot;\xi,\theta):\mathcal{I}\rightarrow\Omega\subset\mathbb{R}^{n}$ denotes the unique solution of the system of ODEs of System (2) with initial condition $\xi\in\Omega$ and we assume it is globally defined, i.e., $\mathcal{I}=[0,+\infty)$ ; and the output $y_{(\xi,\theta)}(t)$ , $t\in\mathcal{S}\subset\mathcal{I}$ , is described by some known function $h(\cdot,\cdot):\Omega\times\Theta\rightarrow\mathbb{R}^{m}$ .

We aim to answer one question: given the output $y_{(x_{0},\theta_{0})}(\cdot)$ for certain $(x_{0},\theta_{0})\in\Omega\times\Theta$ , can we uniquely determine $x_{0}$ (given $\theta_{0}$ ), $\theta_{0}$ (given $x_{0}$ ) or both of them?

We present now formal definitions of the identifiability and observability concepts mentioned in the Introduction. As previously mentioned, some parameter-dependent initial conditions may not allow to distinguish different parameters. Hence, we will consider, for each $\theta\in\Theta$ , a subset $\Omega_{\theta}\subset\Omega$ positively invariant w.r.t. the ODE system in System (2), and the set $\Gamma_{\Theta}=\bigcup_{\theta\in\Theta}\Omega_{\theta}\times\{\theta\}$ . Notice that, if $\Omega_{\theta}=\Omega$ , for each $\theta\in\Theta$ , then $\Gamma_{\Theta}=\Omega\times\Theta$ .

Definition 1 (Identifiability in a time set).

System (2) is identifiable on $\Theta$ in $\mathcal{S}\subset\mathcal{I}$ with initial conditions consistent with the family $\{\Omega_{\theta}\}_{\theta\in\Theta}$ if, for any $\theta_{0}\in\Theta$ , for any $\xi\in\Omega_{\theta_{0}}$ , a different $\theta_{0}^{\prime}\in\Theta$ produces a different output at some time $t\in\mathcal{S}$ , i.e.,

h(x(t;\xi,\theta_{0}),\theta_{0})\neq h(x(t;\xi,\theta_{0}^{\prime}),\theta_{0}^{\prime}).

Equivalently, if $h(x(t;\xi,\theta_{0}),\theta_{0})=h(x(t;\xi,\theta_{0}^{\prime}),\theta_{0}^{\prime})$ , for all $t\in\mathcal{S}$ and any $\xi\in\Omega_{\theta_{0}}$ , then $\theta_{0}=\theta_{0}^{\prime}$ . Moreover, if $\Omega_{\theta}=\Omega$ , for each $\theta\in\Theta$ , we say that System (2) is identifiable on $\Theta$ in $\mathcal{S}\subset\mathcal{I}$ with initial conditions in $\Omega$ .

A typical initial condition which may not allow identifiability are the equilibrium points that depend on the parameters, in the case that different parameter vectors generate the same constant observation. This phenomenon will be illustrated in Section 4.

Definition 2 (Observability in a time set).

System (2) is observable on $\Omega$ in $\mathcal{S}\subset\mathcal{I}$ with parameters in $\Theta$ whether, for any $\theta\in\Theta$ , any two different $x_{0},x_{0}^{\prime}\in\Omega$ produce a different output at some time $t\in\mathcal{S}$ , i.e.,

h(x(t;x_{0},\theta),\theta)\neq h(x(t;x_{0}^{\prime},\theta),\theta).

Equivalently, if $h(x(t;x_{0},\theta),\theta)=h(x(t;x_{0}^{\prime},\theta),\theta)$ , for all $t\in\mathcal{S}$ and any $\theta\in\Theta$ , then $x_{0}=x_{0}^{\prime}$ .

Definition 3 (Identifiability).

System (2) is identifiable on $\Theta$ with initial conditions consistent with the family $\{\Omega_{\theta}\}_{\theta\in\Theta}$ if Definition 1 is satisfied for $\mathcal{S}=\mathcal{I}$ . Moreover, if $\Omega_{\theta}=\Omega$ , for each $\theta\in\Theta$ , we say that System (2) is identifiable on $\Theta$ with initial conditions in $\Omega$ .

Definition 4 (Observability).

System (2) is observable on $\Omega$ with parameters in $\Theta$ if Definition 2 is satisfied for $\mathcal{S}=\mathcal{I}$ .

Notice that identifiability (resp. observability) in a time set $\mathcal{S}$ implies identifiability (resp. observability) in any time set $\mathcal{A}$ such that $\mathcal{S}\subset\mathcal{A}$ . However, the converse is not necessarily true, as demonstrated by the following examples.

Example 1. Consider the following system:

\left\{\begin{array}[]{ccl}x^{\prime}(t)&=&k,\quad x(0)=\xi,\\ y_{(\xi,k)}(t)&=&\max(1,x(t)),\end{array}\right.

(3)

where $k\in\Theta=(0,1]$ is unknown. Let $\Omega=\Omega_{k}=[0,\infty)$ , $k\in\Theta$ , be the positively invariant set w.r.t. the ODE given in (3), $\xi\in\Omega$ , $\mathcal{I}=[0,\infty)$ , and $f(x,k)=k$ (which is Lipschitz-continuous w.r.t. $x$ and continuous w.r.t. $k$ ). The unique solution of the ODE is $x(t)=kt+\xi$ , $t\geq 0.$ The system is identifiable on $\Theta$ in $\mathcal{I}$ with initial conditions in $\Omega$ , since, for any $\xi\in\Omega$ , considering different $k_{1},k_{2}\in\Theta$ , we have that

\max(1,k_{1}t+\xi)=k_{1}t+\xi\neq k_{2}t+\xi=\max(1,k_{2}t+\xi),

for all $t\geq 1/\min\{k_{1},k_{2}\}$ . However, if $\mathcal{S}=[0,1/2]\subset\mathcal{I}$ , $\xi\in[0,1/2]$ and $k_{1},k_{2}\in\Theta$ , with $k_{1}\neq k_{2}$ , then

\max(1,k_{1}t+\xi)=1=\max(1,k_{2}t+\xi),

for any $t\in\mathcal{S}$ , and hence the system is not identifiable on $\Theta$ in $\mathcal{S}$ with initial conditions in $\Omega$ . ∎

Example 2. We consider the same case as the one shown in Example 4. The system is observable on $\Omega$ in $\mathcal{I}$ with parameters in $\Theta$ , since, for any $k\in\Theta$ , considering different $x_{0},x_{0}^{\prime}\in\Omega$ , we have that

\max(1,kt+x_{0})=kt+x_{0}\neq kt+x_{0}^{\prime}=\max(1,kt+x_{0}^{\prime}),

for all $t\geq 1/k$ . However, if $\mathcal{S}=[0,1/2]\subset\mathcal{I}$ and $x_{0},x_{0}^{\prime}\in[0,1/2]\subset\Omega$ , with $x_{0}\neq x_{0}^{\prime}$ , then

\max(1,kt+x_{0})=1=\max(1,kt+x_{0}^{\prime}),

for any $k\in\Theta$ , $t\in\mathcal{S}$ , and hence the system is not observable on $\Omega$ in $\mathcal{S}$ with parameters in $\Theta$ . ∎

A natural extension of this framework consists in treating the parameters as part of the states of an augmented system. If one extends the dynamics with $\dot{\theta}=0$ , then both the identifiability and observability properties can be studied as a particular case of observability in higher dimension: If the extended system is observable, then the original system is both observable and identifiable. However, the reverse implication is not generally true. This is related to the joined observability-identifiability of a system (see [13, 63]).

We are now going to consider that we need to recover both $x_{0}$ and $\theta_{0}$ .

Definition 5 (Joined observability-identifiability in a time set).

System (2) is jointly observable-identifiable on $\Gamma_{\Theta}$ in $\mathcal{S}\subset\mathcal{I}$ if, for any $(x_{0},\theta_{0})\in\Gamma_{\Theta}$ , a different $(x_{0}^{\prime},\theta_{0}^{\prime})\in\Omega\times\Theta$ produces a different output at some time $t\in\mathcal{S}$ , i.e.,

h(x(t;x_{0},\theta_{0}),\theta_{0})\neq h(x(t;x_{0}^{\prime},\theta_{0}^{\prime}),\theta_{0}^{\prime}).

Equivalently, if $h(x(t;x_{0},\theta_{0}),\theta_{0})=h(x(t;x_{0}^{\prime},\theta_{0}^{\prime}),\theta_{0}^{\prime})$ for all $t\in\mathcal{S}$ , this implies that $x_{0}=x_{0}^{\prime}$ and $\theta_{0}=\theta_{0}^{\prime}$ .

Definition 6 (Joined observability-identifiability).

System (2) is jointly observable-identifiable on $\Gamma_{\Theta}$ if Definition 5 is satisfied for $\mathcal{S}=\mathcal{I}$ .

Furthermore, joined observability-identifiability in a time set $\mathcal{S}$ implies joined observability-identifiability in any time set $\mathcal{A}$ such that $\mathcal{S}\subset\mathcal{A}$ .

Note that recovering the parameters and initial condition from the data up to time $t$ is equivalent to reconstructing the parameters and state at current time $t$ , because the dynamics is deterministic and reversible.

Our main focus will be on joined observability-identifiability; however, we will present the results in a way that they may be applicable separately to observability or identifiability. In particular, as commented before, it is clear that joined observability-identifiability implies both observability and identifiability independently.

In the following, we recall the well-known concept of Lie derivative, which is commonly used for studying observability and identifiability. We briefly review it here for completeness.

3.2 About Lie derivatives

Let $\mathbb{N}_{0}=\mathbb{N}\cup\{0\}$ . Consider $\varphi\in\mathcal{C}^{d}(\Omega;\mathbb{R}^{m}),\ F\in\mathcal{C}^{\max\{0,d-1\}}(\Omega;\mathbb{R}^{n})$ , for some $d\in\mathbb{N}\cup\{0\}$ , and the map

\begin{array}[]{lccl}\mathcal{L}_{F,\varphi,d}:&\Omega&\longrightarrow&\mathbb{R}^{m\times(d+1)}\\ &\xi&\longmapsto&\left(L_{F}^{0}\varphi(\xi),L_{F}^{1}\varphi(\xi),\dots,L^{d}_{F}\varphi(\xi)\right),\end{array}

where $L_{F}^{0}\varphi=\varphi$ and, if $d\geq 1$ , $L^{1}_{F}\varphi\in\mathcal{C}^{d-1}(\Omega;\ \mathbb{R}^{m})$ is the Lie derivative of $\varphi$ with respect to the vector field $F$ , i.e., for any $\xi\in\Omega$ ,

L_{F}^{1}\varphi(\xi)=\mathrm{D}\varphi(\xi)F(\xi),

where $\mathrm{D}\varphi$ is the $m\times n$ Jacobian matrix of $\varphi$ . In particular, if $z(\cdot;\xi)$ is a solution of

\dot{z}(t;\xi)=F(z(t;\xi)),\quad z(0;\xi)=\xi,

(4)

such that $\Omega$ is positively invariant with respect to this system, then, for any $t\geq 0$ ,

L_{F}^{1}\varphi(z(t;\xi))=\mathrm{D}\varphi(z(t;\xi))F(z(t;\xi))=\dfrac{\mathrm{d}}{\mathrm{d}t}\varphi(z(t;\xi)).

(5)

Then,

L_{F}^{1}\varphi(\xi)=\left.\dfrac{\mathrm{d}}{\mathrm{d}t}\varphi(z(t;\xi))\right|_{t=0}.

(6)

If $d\geq 2$ , we also define $L^{k}_{F}\varphi\in\mathcal{C}^{d-k}(\Omega;\mathbb{R}^{m})$ , for $k\in\{2,\dots,d\}$ , by recursion as follows:

L^{k}_{F}\varphi=L_{F}^{1}\left(L^{k-1}_{F}\varphi\right).

Lemma 1.

Let $\varphi\in\mathcal{C}^{d}(\Omega;\mathbb{R}^{m})$ , $d\in\mathbb{N}_{0}$ , $\xi\in\Omega$ , and consider System (4), where $F\in\mathcal{C}^{\max\{0,d-1\}}(\Omega;\mathbb{R}^{n})$ and $\Omega$ is positively invariant with respect to the system. Then, for $t\geq 0$ , $k\in\{0,\dots,d\}$ ,

L^{k}_{F}\varphi(z(t;\xi))=\dfrac{\mathrm{d}^{k}}{\mathrm{d}t^{k}}\varphi(z(t;\xi))\ \text{ and, in particular, }\ L^{k}_{F}\varphi(\xi)=\left.\dfrac{\mathrm{d}^{k}}{\mathrm{d}t^{k}}\varphi(z(t;\xi))\right|_{t=0}.

Proof.

For any $\xi\in\Omega$ , $t\geq 0$ , if $d\geq 1$ , we know that (5)-(6) is true. By induction, if $d\geq 2$ , assume that

L^{k-1}_{F}\varphi(z(t;\xi))=\dfrac{\mathrm{d}^{k-1}}{\mathrm{d}t^{k-1}}\varphi(z(t;\xi))\text{ is true for some $k\in\{2,\dots,d\}$, and hence }L^{k-1}_{F}\varphi(\xi)=\left.\dfrac{\mathrm{d}^{k-1}}{\mathrm{d}t^{k-1}}\varphi(z(t;\xi))\right|_{t=0}.

Then, for $k$ , we have

\dfrac{\mathrm{d}^{k}}{\mathrm{d}t^{k}}\varphi(z(t;\xi))=\dfrac{\mathrm{d}}{\mathrm{d}t}L_{F}^{k-1}\varphi(z(t;\xi))=\mathrm{D}L_{F}^{k-1}\varphi(z(t;\xi))F(z(t;\xi))=L_{F}^{k}\varphi(z(t;\xi)),\text{ and }L_{F}^{k}\varphi(\xi)=\left.\dfrac{\mathrm{d}^{k}}{\mathrm{d}t^{k}}\varphi(z(t;\xi))\right|_{t=0},

as we wanted to prove. ∎

Let us now, given $\theta\in\Theta$ , denote $h_{\theta}(\cdot)=h(\cdot,\theta)$ and $f_{\theta}(\cdot)=f(\cdot,\theta)$ , and assume $h_{\theta}\in\mathcal{C}^{d}(\Omega;\mathbb{R}^{m}),\ f_{\theta}\in\mathcal{C}^{\max\{0,d-1\}}(\Omega;\mathbb{R}^{n})$ , for some $d\in\mathbb{N}\cup\{0\}$ . Then, $y_{(\xi,\theta)}\in\mathcal{C}^{d}(\mathcal{I};\mathbb{R}^{m})$ . Indeed, for all $t\geq 0$ ,

\dot{y}_{(\xi,\theta)}(t)=\dfrac{\mathrm{d}}{\mathrm{d}t}h(x(t;\xi,\theta))=\mathrm{D}h_{\theta}(x(t;\xi,\theta))f_{\theta}(x(t;\xi,\theta))\implies\dot{y}_{(\xi,\theta)}(t)=L^{1}_{f_{\theta}}h_{\theta}(x(t;\xi,\theta)).

By Lemma 1, we have, for all $t\geq 0$ and $k\in\{1,\dots,d\}$ ,

y^{(k)}_{(\xi,\theta)}(t)=\dfrac{\mathrm{d}^{k}}{\mathrm{d}t^{k}}h(x(t;\xi,\theta))=L_{f_{\theta}}^{k}h_{\theta}(x(t;\xi,\theta))\text{ and, hence, }y^{(k)}_{(\xi,\theta)}(0)=\left.\dfrac{\mathrm{d}^{k}}{\mathrm{d}t^{k}}h(x(t;\xi,\theta))\right|_{t=0}=L_{f_{\theta}}^{k}h_{\theta}(\xi),

where we denote $y^{(k)}=\dfrac{\mathrm{d}^{k}y}{\mathrm{d}t^{k}}$ , $k\in\mathbb{N},$ when $y\in\mathcal{C}^{k}(\mathcal{I})$ . Then, denoting $y^{(0)}=y$ , define $\mathcal{L}_{f_{\theta},h_{\theta},d}:\Omega\rightarrow\mathbb{R}^{m\times(d+1)}$ as

\mathcal{L}_{f_{\theta},h_{\theta},d}(\xi)=\left(y^{(0)}_{(\xi,\theta)}(0),\dots,y^{(d)}_{(\xi,\theta)}(0)\right).

Roughly speaking, the Lie derivatives of an output of the considered system allow to express the time derivatives of the output not as functions of time, but as functions of the state vector of the system.

In the following, we denote $y_{(\xi,\theta)}$ , $\dot{y}_{(\xi,\theta)}$ and $y^{(k)}_{(\xi,\theta)}$ as $y$ , $\dot{y}$ and $y^{(k)}$ , respectively, when the context is clear.

Let us now present our main results along with the algorithms to recover both the parameters and the initial condition.

3.3 Main results

We revisit and extend some results from [44] to establish sufficient hypotheses for System (2) to be observable, identifiable or jointly observable-identifiable. These results generalize classical approaches (see [44] for a thorough review).

We start recalling a classical result on observability based on Lie derivatives ([13, 24, 31, 47]). This result is of particular interest in the nonlinear context (for which the Cayley-Hamilton theorem is not available, see [13, Section 2.2]). We present the result and a short proof adapted to our framework.

Theorem 1.

Let $\Omega\subset\mathbb{R}^{n}$ positively invariant with respect to the ODE system given in (2), $h_{\theta,i}\in\mathcal{C}^{d_{i}}(\Omega;\mathbb{R}^{m})$ , for some $d_{i}\in\mathbb{N}\cup\{0\}$ , $i\in\{1,\dots,m\}$ , with $m$ the number of scalar outputs, and $f_{\theta}\in\mathcal{C}^{d-1}(\Omega;\mathbb{R}^{n})$ , $d=\max\{1,d_{1},\dots,d_{m}\}$ , for any $\theta\in\Theta$ . If

\begin{array}[]{lclcl}\mathcal{L}_{f_{\theta},h_{\theta},\{d_{1},\dots,d_{m}\}}&:&\Omega&\rightarrow&\mathbb{R}^{m+d_{1}+\dots+d_{m}}\\ &&\xi&\mapsto&\left(\mathcal{L}_{f_{\theta},h_{\theta,1},d_{1}}(\xi),\dots,\mathcal{L}_{f_{\theta},h_{\theta,m},d_{m}}(\xi)\right)\end{array}

is injective in $\Omega$ , then System (2) is observable on $\Omega$ in any semi-open interval $[a,b)\subset\mathcal{I}$ with parameters in $\Theta$ .

Proof.

Given $\theta\in\Theta$ , let $x_{0},x_{0}^{\prime}\in\Omega$ such that

h(x(t;x_{0},\theta),\theta)=h(x(t;x_{0}^{\prime},\theta),\theta),\quad\forall\,t\in[a,b).

Given $\tilde{t}\in[a,b)$ , let $\xi=x(\tilde{t};x_{0},\theta)$ and $\xi^{\prime}=x(\tilde{t};x_{0}^{\prime},\theta).$ Since $\Omega$ is positively invariant w.r.t. the ODE system given in (2), $\xi$ , $\xi^{\prime}\in\Omega$ . Then,

h(x(t;\xi,\theta),\theta)=h(x(t;\xi^{\prime},\theta),\theta),\quad\forall\,t\in[0,b-\tilde{t}).

Notice that $0\in[a-\tilde{t},b-\tilde{t})$ and $[0,b-\tilde{t})\subset\mathcal{I}$ is non-empty. This implies that

y^{(k)}_{(\xi,\theta),i}(t)=y^{(k)}_{(\xi^{\prime},\theta),i}(t),\quad\forall\,t\in[0,b-\tilde{t}),\ k\in\{0,\dots,d_{i}\},\ i\in\{1,\dots,m\},

where $y^{(k)}_{(\xi,\theta)}=\left(y^{(k)}_{(\xi,\theta),1},\dots,y^{(k)}_{(\xi,\theta),m}\right)$ . In particular, $y_{(\xi,\theta),i}^{(k)}(0)=y_{(\xi^{\prime},\theta),i}^{(k)}(0)$ , $k\in\{0,\dots,d_{i}\}$ , $i\in\{1,\dots,m\}$ . Then,

\mathcal{L}_{f_{\theta},h_{\theta},\{d_{1},\dots,d_{m}\}}(\xi)=\mathcal{L}_{f_{\theta},h_{\theta},\{d_{1},\dots,d_{m}\}}(\xi^{\prime}).

Since $\mathcal{L}_{f_{\theta},h_{\theta},\{d_{1},\dots,d_{m}\}}$ is injective in $\Omega$ , for any $\theta\in\Theta$ , and $\xi,\xi^{\prime}\in\Omega$ due to the positive invariance, this implies that $\xi=\xi^{\prime}$ , i.e., $x(\tilde{t};x_{0},\theta)=x(\tilde{t};x_{0}^{\prime},\theta)$ . Due to the uniqueness of solutions of System (2), this implies that $x_{0}=x_{0}^{\prime}$ .

Hence, System (2) is observable on $\Omega$ in any $[a,b)\subset\mathcal{I}$ with parameters in $\Theta$ . ∎

For $(\xi,\theta)\in\Gamma_{\Theta}$ , $\{d_{i}\}_{i=1}^{m}\subset\mathbb{N}_{0}$ , we define the notation

\mathbf{y}_{(\xi,\theta)}^{(d_{1},\dots,d_{m})}(t)\coloneqq\left(y_{(\xi,\theta),1}^{(0)}(t),\dots,y_{(\xi,\theta),1}^{(d_{1})}(t),\dots,y_{(\xi,\theta),m}^{(0)}(t),\dots,y^{(d_{m})}_{(\xi,\theta),m}(t)\right),

and may use the shorthand notation $\mathbf{y}^{(d_{1},\dots,d_{m})}$ when there is no ambiguity.

The following is the main result in [44], adapted to Definition 1, which provides a sufficient condition to ensure recoverability of parameters.

Theorem 2.

Let $h_{\theta,i}\in\mathcal{C}^{d_{i}^{\prime}}(\Omega;\mathbb{R})$ , for some $d_{i}^{\prime}\in\mathbb{N}_{0}$ ¹¹1We use a different notation with respect to Theorem 1 to explicitly state that these orders may be different., $i\in\{1,\dots,m\}$ , and $f_{\theta}\in\mathcal{C}^{d^{\prime}-1}(\Omega;\mathbb{R}^{n})$ , $d^{\prime}=\max\{1,d_{1}^{\prime},\dots,d_{m}^{\prime}\}$ , for any $\theta\in\Theta$ . Consider $\mathcal{D}\subset\mathbb{R}^{d_{1}^{\prime}+\dots+d_{m}^{\prime}+m}$ such that, for all $(\xi,\theta)\in\Gamma_{\Theta}$ , ${\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}(t)\in\mathcal{D}}$ , for all $t\in\mathcal{I}$ . Let $\mathcal{S}\subset\mathcal{I}$ be such that every connected component of $S$ contains an open interval. Assume there exist $g:\mathcal{D}\rightarrow\mathbb{R}^{q+p}$ and $r:\Theta\rightarrow\mathbb{R}^{q}$ , for some $q,p\in\mathbb{N}$ , satisfying:

(C1)

$g=(g_{1,0},\dots,g_{1,q_{1}},\dots,g_{p,0},\dots,g_{p,q_{p}})$ and $r=(r_{1,1},\dots,r_{1,q_{1}},\dots,r_{p,1},\dots,r_{p,q_{p}})$ , with $q_{1}+\dots+q_{p}=q$ , satisfy that

{g_{j,0}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})=\sum_{l=1}^{q_{j}}r_{j,l}(\theta)g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})},

(7)

in $\mathcal{S}$ , for all $j\in\{1,\dots,p\}$ , for any $(\xi,\theta)\in\Gamma_{\Theta}$ ,

(C2)

$r$ is injective, and
(C3)

for any $j\in\{1,\dots,p\}$ , $(\xi,\theta)\in\Gamma_{\Theta}$ , we have that $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}(t))$ , $l\in\{1,\dots,q_{j}\}$ , are linearly independent functions with respect to $t\in\mathcal{S}$ .

Then, System (2) is identifiable on $\Theta$ in $\mathcal{S}$ with initial conditions consistent with the family $\{\Omega_{\theta}\}_{\theta\in\Theta}$ .

Proof.

Given $\theta_{0}\in\Theta$ , $\xi\in\Omega_{\theta_{0}}$ , let $\theta_{0}^{\prime}\in\Theta$ such that

h_{\theta_{0}}(x(t;\xi,\theta_{0}))=h_{\theta_{0}^{\prime}}(x(t;\xi,\theta_{0}^{\prime})),\quad\forall\,t\in\mathcal{S},

i.e.

y_{(\xi,\theta_{0})}(t)=y_{(\xi,\theta_{0}^{\prime})}(t),\quad\forall\,t\in\mathcal{S}.

Then, since every connected part of $\mathcal{S}$ contains some open interval, this implies that

y^{(k)}_{(\xi,\theta_{0}),i}(t)=y^{(k)}_{(\xi,\theta_{0}^{\prime}),i}(t),\quad\forall\,t\in\mathcal{S},\,k\in\{0,\dots,d_{i}^{\prime}\},\,i\in\{1,\dots,m\}.

Since

g_{j,0}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(\xi,\theta_{0})})-g_{j,0}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(\xi,\theta_{0}^{\prime})})\equiv 0

in $\mathcal{S}$ , we obtain, from (7),

\sum_{l=1}^{q_{j}}\big{(}r_{j,l}(\theta_{0})-r_{j,l}(\theta_{0}^{\prime})\big{)}g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(\xi,\theta_{0})})\equiv 0,

in $\mathcal{S},$ for every $j\in\{1,\dots,p\}$ . Given the linear independence of $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})$ , $l\in\{1,\dots,q_{j}\}$ , in $\mathcal{S}\subset\mathcal{I}$ , for every $j\in\{1,\dots,p\}$ , $(\xi,\theta)\in\Gamma_{\Theta}$ , then $r(\theta_{0})=r(\theta_{0}^{\prime})$ . Since $r$ is injective, we have $\theta_{0}=\theta_{0}^{\prime}$ . Hence, System (2) is identifiable on $\Theta$ in $\mathcal{S}\subset\mathcal{I}$ with initial conditions consistent with the family $\{\Omega_{\theta}\}_{\theta\in\Theta}$ . ∎

Remark 1.

It is important to remark that Theorem 2 gives only sufficient conditions for identifiability, but in Section 4.2 we will see an example of a system which is identifiable and does not satisfy these conditions. Nevertheless, as we will see in Algorithm 3, the conditions presented in Theorems 1 and 2 will allow us to have joined observability-identifiability of System (2). In [15, Proposition 2.1], the authors propose a result similar to Theorem 2, giving conditions for it to be also necessary, in a context of rational equations and under the availability of observations at initial time.

Therefore, given a system of ODEs, along with some (partial) observations, we can check the hypotheses in Theorems 1 and 2 in order to determine the observability, identifiability or joined observability-identifiability of our model and recover the initial condition and parameter vector. Notice that we need to be careful when choosing the positively invariant sets we are going to work with, so that the hypotheses mentioned above are satisfied. Next, we present a theorem which shows that, under certain conditions, given $y_{(x_{0},\theta_{0})}$ in some time set, we can recover $(x_{0},\theta_{0})\in\Gamma_{\Theta}$ .

Theorem 3.

For each $\theta\in\Theta$ , let $\Omega_{\theta}\subset\Omega$ be positively invariant with respect to the ODE system in (2). Let $(x_{0},\theta_{0})\in\Gamma_{\Theta}$ . Assume we know $y_{(x_{0},\theta_{0})}$ in $\mathcal{S}\subset\mathcal{I}$ , $\mathcal{S}$ such that every connected component contains an open interval. If the hypotheses of Theorems 1 and 2 are satisfied, then System (2) is jointly observable-identifiable in $\mathcal{S}$ and we can reconstruct the pair $(x_{0},\theta_{0})$ univocally using the values of $y_{(x_{0},\theta_{0})}$ and its derivatives at, at most, $q+1=q_{1}+\dots+q_{p}+1$ suitable values of $t\in\mathcal{S}$ .

Proof.

Let $\phi_{j,l}(t)=g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(x_{0},\theta_{0})}(t))$ , $t\in\mathcal{S}$ , $l\in\{0,\dots,q_{j}\},\ j\in\{1,\dots,p\}$ . By hypothesis, $\{\phi_{j,l}\}_{l=1}^{q_{j}}$ are linearly independent in $\mathcal{S}$ , for each $j\in\{1,\dots,p\}$ . Then, for every $j\in\{1,\dots,p\}$ , there exist (see [44, Lemma 1]) $q_{j}$ different time instants $\{t_{j,\ell}\}_{\ell=1}^{q_{j}}\subset\mathcal{S}$ such that

\det((\phi_{j,l}(t_{j,\ell}))_{l,\ell=1,\dots,q_{j}})\neq 0.

Then, there exists a unique solution $\sigma$ to

\left(\phi_{j,1}(t_{j,\ell})\ \dots\ \phi_{j,q_{j}}(t_{j,\ell})\right)\sigma=\phi_{j,0}(t_{j,\ell}),\ \forall\ \ell\in\{1,\dots,q_{j}\}.

(8)

Attending to (7), it satisfies $\sigma_{j,l}=r_{j,l}(\theta_{0})$ , $l\in\{1,\dots,q_{j}\},\ j\in\{1,\dots,p\}$ . Since $r$ is injective, such that $r^{-1}:r(\Theta)\rightarrow\Theta$ , we can recover our original parameter vector

\theta_{0}=r^{-1}(\sigma_{1,1},\dots,\sigma_{p,q_{p}}).

Finally, to recover $x_{0}$ , take some time $\tilde{t}\in\mathcal{S}$ , which can be some $\tilde{t}\in\{t_{1,1},\dots,t_{p,q_{p}}\}$ . Since $\mathcal{L}_{f_{\theta_{0}},h_{\theta_{0}},\{d_{1},\dots,d_{m}\}}$ is injective in $\Omega$ , there exists a unique $\tilde{\xi}\in\Omega$ such that

\tilde{\xi}=\mathcal{L}^{-1}_{f_{\theta_{0}},h_{\theta_{0}},\{d_{1},\dots,d_{m}\}}{(\mathbf{y}^{(d_{1},\dots,d_{m})}_{(x_{0},\theta_{0})}(\tilde{t}))},

noticing that $\mathbf{y}^{(d_{1},\dots,d_{m})}_{(x_{0},\theta_{0})}(\tilde{t})=\mathbf{y}^{(d_{1},\dots,d_{m})}_{(\tilde{\xi},\theta_{0})}(0)$ . Since it is unique, it must satisfy $\tilde{\xi}\in\Omega_{\theta_{0}}$ . We can recover $x_{0}\in\Omega_{\theta_{0}}$ integrating backwards the ODE system in System (2) knowing $\tilde{\xi}$ , $\tilde{t}$ and $\theta_{0}$ (since $f_{\theta_{0}}$ is Lipschitz in $\Omega$ and, in particular, in $\Omega_{\theta_{0}}$ positively invariant w.r.t. the ODE system of (2)). If $0\in\mathcal{S}$ , we directly choose $\tilde{t}=0$ .

This is, we have recovered $\theta_{0}$ and $x_{0}$ from the data univocally knowing $y_{(x_{0},\theta_{0})}$ in $\mathcal{S}$ , using its value and the value of its derivatives at $q+1$ (at most) different time instants. ∎

In the following Algorithm 3, we present a procedure to recover the unknowns given some observations, based on the proof of Theorem 3.

Algorithm 1. Assume that we know $y_{(x_{0},\theta_{0})}$ (satisfying System (2)) in $\mathcal{S}\subset\mathcal{I}$ , such that every connected of $\mathcal{S}$ component contains an open interval, and $x_{0}\in\Omega$ and/or $\theta_{0}\in\Theta$ are unknown. The procedure to recover the unknowns is the following:

Step 1.

If $x_{0}$ is unknown, find $\Omega_{1}\subset\Omega$ positively invariant with respect to the ODE system given in (2) such that $x_{0}\in\Omega_{1}$ , and $d_{1},\dots,d_{m}\in\mathbb{N}_{0}$ such that, for any $\theta\in\Theta$ , the following function is injective in $\Omega_{1}$ :

$\mathcal{L}_{\theta,\{d_{1},\dots,d_{m}\}}:\xi\mapsto\mathbf{y}^{(d_{1},\dots,d_{m})}_{(\xi,\theta)}(0).$

Note: According to Theorem 1, this ensures that System (2) is observable on $\Omega_{1}$ in any semi-open interval $[a,b)\subset\mathcal{I}$ with parameters in $\Theta$ .
Step 2.

If $\theta_{0}$ is unknown, find $\Omega_{2,\theta}\subset\Omega$ , for any $\theta\in\Theta$ , positively invariant with respect to the ODE system given in (2) such that $x_{0}\in\Omega_{2,\theta_{0}}$ , and maps $g:\mathcal{D}\rightarrow\mathbb{R}^{q+p}$ and $r:\Theta\rightarrow\mathbb{R}^{p}$ , for some suitable $\mathcal{D}\subset\mathbb{R}^{d_{1}^{\prime}+\dots+d_{m}^{\prime}+m}$ , $d_{1}^{\prime},\dots,d_{m}^{\prime}\in\mathbb{N}_{0}$ , $q,p\in\mathbb{N}$ , such that (C1), (C2) and (C3) of Theorem 2 are satisfied.

Note: According to Theorem 2, this ensures that System (2) is identifiable on $\Theta$ in $\mathcal{S}$ with initial conditions consistent with the family $\{\Omega_{2,\theta}\}_{\theta\in\Theta}$ .
Step 3.

If we have defined $\Omega_{1}$ and $\Omega_{2,\theta}$ , set $\Omega_{\theta}=\Omega_{1}\cap\Omega_{2,\theta}$ , for all $\theta\in\Theta$ . If we have only defined $\Omega_{1}$ , redefine $\Omega=\Omega_{1}$ . If we have only defined $\Omega_{2,\theta}$ , set $\Omega_{\theta}=\Omega_{2,\theta}$ , for every $\theta\in\Theta$ .
Step 4.

If $\theta_{0}$ is unknown, for each $j\in\{1,\dots,p\}$ , find $q_{j}$ different time instants $t_{j,1},\dots,t_{j,q_{j}}\in\mathcal{S}$ such that

$\det\left(\phi_{j,l}(t_{j,\ell})_{l,\ell=1,\dots,q_{j}}\right)\neq 0,$

where $\phi_{j,l}=g_{j,l}\left(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(x_{0},\theta_{0})}\right)$ .

Note: These time instants exist due to condition (C3) of Theorem 2, which is verified in Step 2, and [44, Lemma 1].

Step 5.

If $\theta_{0}$ is unknown, for each $j\in\{1,\dots,p\}$ , solve the following linear system which, according to Step 4, has a unique solution $\sigma_{j}=(\sigma_{j,1},\dots,\sigma_{j,q_{j}})$ :

\left(\begin{array}[]{ccc}\phi_{j,1}(t_{j,1})&\cdots&\phi_{j,q_{j}}(t_{j,1})\\ \vdots&\ddots&\vdots\\ \phi_{j,1}(t_{j,q_{j}})&\cdots&\phi_{j,q_{j}}(t_{j,q_{j}})\end{array}\right)\left(\begin{array}[]{c}\sigma_{j,1}\\ \vdots\\ \sigma_{j,q_{j}}\end{array}\right)=\left(\begin{array}[]{c}\phi_{j,0}(t_{j,1})\\ \vdots\\ \phi_{j,0}(t_{j,q_{j}})\end{array}\right),

(9)

where $\phi_{j,0}=g_{j,0}\left(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(x_{0},\theta_{0})}\right)$ .

Step 6.

If $\theta_{0}$ is unknown, recover $\theta_{0}$ solving the equation

$r(\theta_{0})=(\sigma_{1},\dots,\sigma_{p}).$

Note: We are taking into account that $r$ is injective and $r(\theta_{0})$ is solution of System (9), according to conditions (C1) and (C2) of Theorem 2, which are verified in Step 2.

Step 7.

If $x_{0}$ is unknown, choose some $\tilde{t}\in\mathcal{S}$ (which can be some $\tilde{t}\in\{t_{1,1},\dots,t_{p,q_{p}}\}$ if we performed Step 4) and solve for $\tilde{\xi}$ the equation

\mathcal{L}_{\theta_{0},\{d_{1},\dots,d_{m}\}}(\tilde{\xi})=\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(x_{0},\theta_{0})}(\tilde{t}).

Note: We are taking into account that $\mathcal{L}_{\theta,\{d_{1},\dots,d_{m}\}}$ is injective in $\Omega$ and, in particular, in $\Omega_{\theta}$ , for any $\theta\in\Theta$ , according to Step 1, and $y_{(x_{0},\theta_{0}),i}^{(k)}(\tilde{t})=y_{(\tilde{\xi},\theta_{0}),i}^{(k)}(0)$ , for all $k\in\{0,\dots,d_{i}\}$ , $i\in\{1,\dots,m\}$ .

Step 8.

If $x_{0}$ is unknown and $\tilde{t}\neq 0$ , integrate backwards System (2) from $\tilde{t}$ to $0$ with initial condition $\tilde{\xi}$ , knowing $\theta_{0}$ , in order to recover $x_{0}$ .

Recall that the time instants required in Step 4 may be anywhere in $\mathcal{S}\subset\mathcal{I}$ , and hence may be difficult to find in practice. Nevertheless, in Lemma (2) we give sufficient hypotheses such that System (2) is identifiable in any semi-open interval $[a,b)\subset\mathcal{S}$ ; this will imply that we will be able to choose this set of time instants in any of these semi-open intervals. Then, an analogous procedure to this one will be presented in Algorithm 3.

Remark 2.

Notice that, in order to be able to recover $(x_{0},\theta_{0})$ following the procedure in Step 4 of Algorithm 3, for each set $\{\phi_{j,1},\dots,\phi_{j,q_{j}}\}$ , $j\in\{1,\dots,p\}$ , we need to find $q_{j}$ different suitable time instants. However, some of these time instants may coincide among different sets of linearly independent functions. Thus, the number of different time instants we need to find is between $\tilde{q}=\max\{q_{1},\dots,q_{p}\}$ and $q=q_{1}+\dots+q_{p}$ , along maybe with $\tilde{t}=0$ , which we can use to recover the initial condition if $0\in\mathcal{S}$ and could be one of the other time instants.

Recall, moreover, that we do not necessarily need $y_{(x_{0},\theta_{0})}(t)$ , for all $t\in\mathcal{S}$ , but it would be sufficient having the values of $y_{(x_{0},\theta_{0})}$ and its derivatives at the aforementioned different time instants, where the order of the derivatives that we need are the same as for Algorithm 3.

Note: Although this idea of choosing a number of suitable time instants has already been commented in [59, Remark 3] and [42], no proof nor more detail is given.

Lemma 2.

Assume the hypotheses of Theorem 2 are satisfied for some $\mathcal{S}\subset\mathcal{I}$ and, for any $j\in\{1,\dots,p\}$ , $l\in\{1,\dots,q_{j}\}$ and $(\xi,\theta)\in\Gamma_{\Theta}$ , the functions $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})$ are analytic on $\mathcal{I}$ . Then, System (2) is identifiable on $\Theta$ in any semi-open interval $[a,b)\subset\mathcal{S}$ with initial conditions consistent with the family $\{\Omega_{\theta}\}_{\theta\in\Theta}$ . If $\mathcal{S}$ is connected, it is sufficient for $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})$ , $j\in\{1,\dots,p\}$ , $l\in\{1,\dots,q_{j}\}$ , to be analytic on $\mathcal{S}$ .

Proof.

Given $\theta_{0}\in\Theta$ , $\xi\in\Omega_{\theta_{0}}$ , let $\theta_{0}^{\prime}\in\Theta$ such that

h_{\theta_{0}}(x(t;\xi,\theta_{0}))=h_{\theta_{0}^{\prime}}(x(t;\xi,\theta_{0}^{\prime})),\quad\forall\,t\in[a,b),

i.e.,

y_{(\xi,\theta_{0})}(t)=y_{(\xi,\theta_{0}^{\prime})}(t),\quad\forall\,t\in[a,b).

Then,

y^{(k)}_{(\xi,\theta_{0}),i}(t)=y^{(k)}_{(\xi,\theta_{0}^{\prime}),i}(t),\quad\forall\,t\in[a,b),

$k\in\{0,\dots,d_{i}^{\prime}\},\ i\in\{1,\dots,m\}.$ Since, for any $j\in\{1,\dots,p\}$ , $l\in\{1,\dots,q_{j}\}$ and $(\xi,\theta)\in\Gamma_{\Theta}$ , the functions $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})$ are analytic in $\mathcal{I}$ , then the function

G_{j,\theta}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})\coloneqq\sum_{l=1}^{q_{j}}r_{j,l}(\theta)g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})

is also analytic in $\mathcal{I}$ . Given that $[a,b)\subset\mathcal{S}$ ,

G_{j,\theta}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})=g_{j,0}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}),

in $[a,b)$ , for any $\theta\in\Theta$ . Since

g_{j,0}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(\xi,\theta_{0})})=g_{j,0}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(\xi,\theta_{0}^{\prime})}),

in $[a,b)$ , we have

G_{j,\theta_{0}}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(\xi,\theta_{0})})=G_{j,\theta_{0}^{\prime}}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(\xi,\theta_{0}^{\prime})}),

in $[a,b)$ . This implies that, for every $j\in\{1,\dots,p\}$ ,

R_{j}\coloneqq\sum_{l=1}^{q_{j}}\big{(}r_{j,l}(\theta_{0})-r_{j,l}(\theta_{0}^{\prime})\big{)}g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(\xi,\theta_{0})})\equiv 0

in $[a,b)$ . Due to the analyticity of $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})$ , $l\in\{1,\dots,q_{j}\}$ , in $\mathcal{I}$ , $R_{j}$ is also analytic in $\mathcal{I}$ . If $R_{j}\equiv 0$ in $[a,b)$ , then $R_{j}\equiv 0$ in $\mathcal{I}$ ([58, Theorem 8.5]). Therefore, since $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})$ , $l\in\{1,\dots,q_{j}\}$ , are linearly independent in $\mathcal{S}\subset\mathcal{I}$ , they are also linearly independent in $\mathcal{I}$ , and hence, for $R_{j}$ to be 0 in $\mathcal{I}$ , for all $j\in\{1,\dots,p\}$ , we need $r(\theta_{0})=r(\theta_{0}^{\prime})$ . Since $r$ is injective, this implies $\theta_{0}=\theta_{0}^{\prime}$ .

If $\mathcal{S}$ is connected, $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})$ , $l\in\{1,\dots,q_{j}\}$ , analytic in $\mathcal{S}$ implies $R_{j}$ is analytic in $\mathcal{S}$ , and $R_{j}\equiv 0$ in $[a,b)\subset\mathcal{S}$ implies $R_{j}\equiv 0$ in $\mathcal{S}$ (if $\mathcal{S}$ is not connected, we can only assure $R_{j}\equiv 0$ in the connected component containing $[a,b)$ ). Thus, we conclude analogously using the linear independence of $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})$ , $l\in\{1,\dots,q_{j}\}$ , in $\mathcal{S}$ . Hence, System (2) is identifiable on $\Theta$ in any $[a,b)\subset\mathcal{S}$ with initial conditions consistent with the family $\{\Omega_{\theta}\}_{\theta\in\Theta}$ . ∎

Taking into account Theorem 1 and Lemma 2, we present a lemma that shows that we may recover $x_{0}$ and $\theta_{0}$ knowing $y_{(x_{0},\theta_{0})}$ only in some $[a,b)\subset\mathcal{I}$ .

Lemma 3.

For each $\theta\in\Theta$ , let $\Omega_{\theta}\subset\Omega$ be positively invariant with respect to the ODE system in (2). Let $(x_{0},\theta_{0})\in\Gamma_{\Theta}$ . Assume we know $y_{(x_{0},\theta_{0})}$ in some $[a,b)\subset\mathcal{I}$ . If the hypotheses of Theorem 1 and Lemma 2 are satisfied for some $\mathcal{S}\subset\mathcal{I}$ such that $[a,b)\subset\mathcal{S}$ , then System (2) is jointly observable-identifiable in $[a,b)$ , and we can reconstruct $(x_{0},\theta_{0})$ univocally using the values of $y_{(x_{0},\theta_{0})}$ and its derivatives at a finite amount of suitable values in $[a,b)$ . In particular, the required number of time points is between $\tilde{q}=\max\{q_{1},\dots,q_{p}\}$ and $q=q_{1}+\dots+q_{p}$ .

Proof.

We recover $\theta_{0}$ analogously to the proof of Theorem 3. We only need to see that, given $j\in\{1,\dots,p\}$ , since the functions $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})$ , $l\in\{1,\dots,q_{j}\}$ , are linearly independent in some $\mathcal{S}\subset\mathcal{I}$ and analytic in $\mathcal{I}$ , they are linearly independent in $[a,b)\subset\mathcal{S}$ . Given $j\in\{1,\dots,p\}$ , assume that $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})$ , $l\in\{1,\dots,q_{j}\}$ , are linearly dependent in $[a,b)$ , i.e., there exist $\{a_{j,l}\}_{l=1}^{q_{j}}\subset\mathbb{R}$ not all of them null such that,

G_{j}(t)=\sum_{l=1}^{q_{j}}a_{j,l}g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}(t))=0,\ \forall\ t\in[a,b).

Since $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})$ , $l\in\{1,\dots,q_{j}\}$ , are analytic in $\mathcal{I}$ , then so it is $G_{j}$ . This implies that (see [58, Theorem 8.5]) $G_{j}\equiv 0$ in $\mathcal{I}$ and, thus, $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})$ , $l\in\{1,\dots,q_{j}\}$ , are linearly dependent in $\mathcal{I}$ , and hence in $\mathcal{S}$ , which is a contradiction.

Moreover, if $\mathcal{S}$ is connected, it is enough asking for $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})$ , $l\in\{1:q_{j}\}$ , analytic in $\mathcal{S}$ , since, then, so it is $G_{j}$ . Therefore, $G_{j}\equiv 0$ in $[a,b)\subset\mathcal{S}$ connected implies (see [58, Theorem 8.5]) that $G_{j}\equiv 0$ in $\mathcal{S}$ , which leads to the same contradiction. Therefore, for each $j\in\{1,\dots,p\}$ , the functions $\{\phi_{j,l}\}_{l=1}^{q_{j}}$ , with $\phi_{j,l}=g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(x_{0},\theta_{0})})$ , $l\in\{1,\dots,q_{j}\}$ , are linearly independent in $[a,b)$ and we can conclude analogously to Theorem 3, along with Remark 2, choosing $t_{1,1},$ $\dots,$ $t_{p,q_{p}}\in[a,b)$ .

On the other hand, to recover the initial condition, we proceed as in the proof of Theorem 3.

Hence, we are able to recover $(x_{0},\theta_{0})$ univocally when knowing $y_{(x_{0},\theta_{0})}$ in $[a,b)$ , using its values and the values of its derivatives at some finite set of time instants in $[a,b)$ ; concretely, between $\tilde{q}$ and $q$ ddistinct suitable time points. ∎

We present in Algorithm 3 a method to recover $x_{0}$ and/or $\theta_{0}$ knowing $y_{(x_{0},\theta_{0})}$ only in some $[a,b)\subset\mathcal{I}$ , based on the proof of previous Lemma 3. In particular, we will need to find between $\tilde{q}=\max\{q_{1},\dots,q_{p}\}$ and $q=q_{1}+\dots+q_{p}$ different time instants in $[a,b)$ when some analyticity properties are satisfied.

Algorithm 2. In a way similar to that of Algorithm 3, we describe here a slightly different constructive method to recover $x_{0}\in\Omega$ and/or $\theta_{0}\in\Theta$ assuming that we know $y_{(x_{0},\theta_{0})}$ (satisfying System (2)) in some $[a,b)\subset\mathcal{I}$ and some analyticity hypotheses are fulfilled. Since the procedure is similar to the one in Algorithm 3, we only state here the steps that change:

Step 2.

If $\theta_{0}$ is unknown, find $\Omega_{2,\theta}\subset\Omega$ , for any $\theta\in\Theta$ , positively invariant with respect to the ODE system given in (2) such that $x_{0}\in\Omega_{2,\theta_{0}}$ , and maps $g:\mathcal{D}\rightarrow\mathbb{R}^{q+p}$ and $r:\Theta\rightarrow\mathbb{R}^{p}$ , for some suitable $\mathcal{D}\subset\mathbb{R}^{d_{1}^{\prime}+\dots+d_{m}^{\prime}+m}$ , $d_{1}^{\prime},\dots,d_{m}^{\prime}\in\mathbb{N}\cup\{0\}$ , $q,p\in\mathbb{N}$ , such that (C1), (C2) and (C3) of Theorem 2 are satisfied, being $\mathcal{S}\subset\mathcal{I}$ in (C2) and (C3) a connected set satisfying $[a,b)\subset\mathcal{S}$ , and such that functions $g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})})$ are analytic in $\mathcal{S}$ , for any $j\in\{1,\dots,p\}$ , $l\in\{1,\dots,q_{j}\}$ , $(\xi,\theta)\in\Gamma_{2,\Theta}=\bigcup_{\theta\in\Theta}\Omega_{2,\theta}\times\{\theta\}$ .

Note: According to Lemma 2, this ensures that System (2) is identifiable on $\Theta$ in any semi-open interval $[a,b)\subset\mathcal{S}$ with initial conditions consistent with the family $\{\Omega_{2,\theta}\}_{\theta\in\Theta}$ .
Step 4.

Same as in Algorithm 3 but such that, for each $j\in\{1,\dots,p\}$ , $t_{j,1},\dots,t_{j,q_{j}}\in[a,b)$ .

Note: These time instants exist due to condition (C3) of Theorem 2 (which is verified in Step 2), the analyticity condition required also in Step 2, and [44, Lemma 1].
Step 7.

Same as in Algorithm 3, but such that $\tilde{t}\in[a,b)$ .

Notice again that we need to find in $[a,b)$ the sets of time instants required in Step 4, and hence again may be difficult to find in practice. However, we can reduce this quantity of time instants to $1$ time if we use higher derivatives of $y$ , as will be shown next.

Remark 3.

Recall that, as in Remark 2, we may not need to know $y_{(x_{0},\theta_{0})}$ in some interval $[a,b)$ , but only the values of this function and its derivatives at the time instants indicated in Step 4 and Step 7 in Algorithm 3, where the needed order of the derivatives is given by the same algorithm.

3.4 Using higher-order derivatives

In this section, we present a methodology which uses higher derivatives of $y$ and may be more convenient in some cases. Some classical differential algebra methodologies use higher-order derivatives and, once they obtain the expressions of these derivatives, they study the identifiability of the system from a differential algebra viewpoint (recall, for example, [26, 28, 52]). In the methodology presented in Lemma 4 and Algorithm 4 we still determine the identifiability of System (2) studying the linear independence, which in several cases will be easier than what is performed in the other methodologies. The use we make of higher-order derivatives allows to reduce to one the number of time instants that we need to find in $[a,b)\subset\mathcal{I}$ to recover the unknown parameters, once we know the system is identifiable. Therefore, notice that using these higher-order derivatives is an alternative to recover the unknowns using less time instants but is not necessary for proving the identifiability of the system, as it is in the aforementioned works.

Lemma 4.

For each $\theta\in\Theta$ , let $\Omega_{\theta}\subset\Omega$ be positively invariant with respect to the ODE system in (2). Let $(x_{0},\theta_{0})\in\Gamma_{\Theta}$ . Assume we know $y_{(x_{0},\theta_{0})}$ in some semi-open interval $[a,b)\subset\mathcal{I}$ , and the hypotheses of Theorem 1 and Lemma 2 are satisfied for some $\mathcal{S}\subset\mathcal{I}$ such that $[a,b)\subset\mathcal{S}$ . If we further assume that $h_{\theta,i}\in\mathcal{C}^{d_{i}^{\prime}+\tilde{q}-1}(\Omega;\mathbb{R})$ , $i\in\{1,\dots,m\}$ , $f_{\theta}\in\mathcal{C}^{d^{\prime}+\tilde{q}-2}(\Omega;\mathbb{R}^{n})$ , for any $\theta\in\Theta$ , being $\tilde{q}=\max\{q_{1},\dots,q_{p}\}$ , then System (2) is jointly observable-identifiable and we can reconstruct $(x_{0},\theta_{0})$ univocally using the value of $y_{(x_{0},\theta_{0})}$ and its derivatives at only one time in $[a,b)$ ; in particular, this time can be almost any $t\in[a,b)$ .

Proof.

Let us take $\phi_{j,l}=g_{j,l}(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(x_{0},\theta_{0})})$ in $[a,b)$ , for all $l\in\{0,\dots,q_{j}\},\ j\in\{1,\dots,p\}$ . Since

\phi_{j,0}(t)=\sum_{l=1}^{q_{j}}r_{j,l}(\theta_{0})\phi_{j,l}(t),\quad\forall\,t\in[a,b),

where $r(\theta_{0})$ does not depend on time, we differentiate each of these equations, for $j\in\{1,\dots,p\}$ , $q_{j}-1$ times such that

\dfrac{\mathrm{d}^{k}\phi_{j,0}}{\mathrm{d}t^{k}}(t)=\sum_{l=1}^{q_{j}}r_{j,l}(\theta_{0})\dfrac{\mathrm{d}^{k}\phi_{j,l}}{\mathrm{d}t^{k}}(t),\quad\forall\ t\in[a,b),\ k\in\{0,\dots,q_{j}-1\},

By the proof of [44, Lemma 3], we know that $\phi_{j,1},\dots,\phi_{j,q_{j}}$ are linearly independent in $[a,b)$ , and hence there exists some $t_{j}\in[a,b)$ (see [5]) such that

W_{j}(t_{j})=\det\left(\left(\dfrac{\mathrm{d}^{k}\phi_{j,l}}{\mathrm{d}t^{k}}(t_{j})\right)_{\begin{subarray}{l}l=1,\dots,q_{j}\\ k=0,\dots,q_{j}-1\end{subarray}}\right)\neq 0,

where $W_{j}$ is the Wrońskian of $\phi_{j,1},\dots,\phi_{j,q_{j}}$ , i.e., there exists a unique solution $\sigma_{j}$ to

\left(\begin{array}[]{ccc}\phi_{j,1}(t_{j})&\cdots&\phi_{j,q_{j}}(t_{j})\\[5.0pt] \dfrac{\mathrm{d}\phi_{j,1}}{\mathrm{d}t}(t_{j})&\cdots&\dfrac{\mathrm{d}\phi_{j,q_{j}}}{\mathrm{d}t}(t_{j})\\[5.0pt] \vdots&\ddots&\vdots\\[5.0pt] \dfrac{\mathrm{d}^{q_{j}-1}\phi_{j,1}}{\mathrm{d}t^{q_{j}-1}}(t_{j})&\cdots&\dfrac{\mathrm{d}^{q_{j}-1}\phi_{j,q_{j}}}{\mathrm{d}t^{q_{j}-1}}(t_{j})\end{array}\right)\left(\begin{array}[]{c}\sigma_{j,1}\\ \sigma_{j,2}\\ \vdots\\ \sigma_{j,q_{j}}\end{array}\right)=\left(\begin{array}[]{c}\phi_{j,0}(t_{j})\\[5.0pt] \dfrac{\mathrm{d}\phi_{j,0}}{\mathrm{d}t}(t_{j})\\[5.0pt] \vdots\\[5.0pt] \dfrac{\mathrm{d}^{q_{j}-1}\phi_{j,0}}{\mathrm{d}t^{q_{j}-1}}(t_{j})\end{array}\right).

Moreover, $W_{j}$ is also analytic, and therefore it can only vanish at isolated points [58, Theorem 8.5], i.e., the previous equation is fulfilled for almost every $t_{j}\in[a,b)$ . Furthermore, since each $W_{j}$ , $j\in\{1,\dots,p\}$ , only vanishes at a countable set of times, then all of them are non-null simultaneously for almost every $t\in[a,b)$ . Hence, we can choose almost any time $\tilde{t}\in[a,b)$ such that the previous system has a unique solution considering $t_{j}=\tilde{t}$ , for every $j\in\{1,\dots,p\}$ .

In a way similar to that of the proof of [44, Theorem 3], we may recover our original parameter vector $\theta_{0}$ as

\theta_{0}=r^{-1}(\sigma_{1,1},\dots,\sigma_{p,q_{p}}).

Having recovered $\theta_{0}$ , we can recover $x_{0}$ in the same way as exposed in the proof of [44, Lemma 3].

Hence, we are able to recover $(x_{0},\theta_{0})$ univocally when knowing $y_{(x_{0},\theta_{0})}(t)$ , $t\in[a,b)$ , using its value and the value of its derivatives at one time in $[a,b)$ , which can be almost any $t\in[a,b)$ . ∎

Therefore, if we ask for higher regularity of the components of $y_{(\xi,\theta)}$ , $(\xi,\theta)\in\Gamma_{\Theta}$ , than required in [44, Lemma 3], we can recover $x_{0}$ and $\theta_{0}$ choosing only one time in $[a,b)$ , instead of a number of distinct time instants between $\tilde{q}=\max\{q_{1},\dots,q_{p}\}\geq 1$ and $q=q_{1}+\dots+q_{p}\geq p$ (see Remark 2). We present next the associated algorithm.

Algorithm 3. In a way similar to those of Algorithm 3 and Algorithm 3, we describe here a different constructive method to recover $x_{0}\in\Omega$ and/or $\theta_{0}\in\Theta$ assuming that we know $y_{(x_{0},\theta_{0})}$ (satisfying System (2)) in some $[a,b)\subset\mathcal{I}$ , some analyticity properties are fulfilled and $y_{(x_{0},\theta_{0})}$ is smooth enough. This method is presented in the proof of Lemma 4, and it requires to choose only one time in $[a,b)$ . However, it needs the use of higher derivatives of $y_{(x_{0},\theta_{0})}$ , which are not necessary in the constructive methods explained in Algorithm 3 and Algorithm 3. We only state here the steps that are different from those in Algorithms 3 and 3:

Step 2.

Same as in Step 2 of Algorithm 3, along with $y_{(\xi,\theta),i}\in\mathcal{C}^{d_{i}^{\prime}+\tilde{q}-1}$ , for every $(\xi,\theta)\in\Gamma_{2,\Theta}$ and each $i\in\{1,\dots,m\}$ , with $\tilde{q}=\max\{q_{1},\dots,q_{p}\}$ .
Step 4.

If $\theta_{0}$ is unknown, for each $j\in\{1,\dots,p\}$ , differentiate with respect to time $q_{j}-1$ times the equation

$\phi_{j,0}=\sum_{l=1}^{q_{j}}r_{j,l}(\theta_{0})\phi_{j,l},$

where $\phi_{j,\ell}=g_{j,\ell}\left(\mathbf{y}^{(d_{1}^{\prime},\dots,d_{m}^{\prime})}_{(x_{0},\theta_{0})}\right)$ , $\ell\in\{0,l\}$ .

Note: We can differentiate due to the analyticity condition on functions $g_{j,l}$ and the regularity conditions on $y$ required in Step 2.

Step 5.

If $\theta_{0}$ is unknown, choose a time $\tilde{t}\in[a,b)$ such that

W_{j}(\tilde{t})=\det\left(\left(\dfrac{\mathrm{d}^{k}\phi_{j,l}}{\mathrm{d}t^{k}}(\tilde{t})\right)_{\begin{subarray}{l}l=1,\dots,q_{j}\\ k=0,\dots,q_{j}-1\end{subarray}}\right)\neq 0,\quad\forall\,j\in\{1,\dots,p\},

(10)

where $W_{j}$ is the Wrońskian of $\{\phi_{j,1},\dots,\phi_{j,q_{j}}\}$ .

Note: This time can be almost any time in $[a,b)$ due to condition (C3) of Theorem 2, which is verified in Step 2, and the analyticity condition required in Step 2.

Step 6.

If $\theta_{0}$ is unknown, for each $j\in\{1,\dots,p\}$ , solve the following linear system which, due to Step 5, has a unique solution $\sigma_{j}=(\sigma_{j,1},\dots,\sigma_{j,q_{j}})$ :

\left(\begin{array}[]{ccc}\phi_{j,1}(\tilde{t})&\cdots&\phi_{j,q_{j}}(\tilde{t})\\[5.0pt] \dfrac{\mathrm{d}\phi_{j,1}}{\mathrm{d}t}(\tilde{t})&\cdots&\dfrac{\mathrm{d}\phi_{j,q_{j}}}{\mathrm{d}t}(\tilde{t})\\[5.0pt] \vdots&\ddots&\vdots\\[5.0pt] \dfrac{\mathrm{d}^{q_{j}-1}\phi_{j,1}}{\mathrm{d}t^{q_{j}-1}}(\tilde{t})&\cdots&\dfrac{\mathrm{d}^{q_{j}-1}\phi_{j,q_{j}}}{\mathrm{d}t^{q_{j}-1}}(\tilde{t})\end{array}\right)\left(\begin{array}[]{c}\sigma_{j,1}\\ \sigma_{j,2}\\ \vdots\\ \sigma_{j,q_{j}}\end{array}\right)=\left(\begin{array}[]{c}\phi_{j,0}(\tilde{t})\\[5.0pt] \dfrac{\mathrm{d}\phi_{j,0}}{\mathrm{d}t}(\tilde{t})\\[5.0pt] \vdots\\[5.0pt] \dfrac{\mathrm{d}^{q_{j}-1}\phi_{j,0}}{\mathrm{d}t^{q_{j}-1}}(\tilde{t})\end{array}\right).

(11)

Step 7.

Same as in Step 6 of Algorithm 3, but considering that $r(\theta_{0})$ is solution of System (11).
Step 8.

Same as in Step 7 of Algorithm 3, but now the time we choose can be the same as $\tilde{t}$ in Steps 5 and 6 of this algorithm.
Step 9.

Same as in Step 8 of Algorithm 3.

Remark 4.

Notice that, as in Remark 2 and Remark 3, we may not need to know $y_{(x_{0},\theta_{0})}$ in some interval $[a,b)$ , but only its value and the values of its derivatives (each component $y_{(x_{0},\theta_{0}),i}$ up to the $\max\{d_{i},d_{i}^{\prime}+\tilde{q}-1\}-$ th derivative) at one time such that Step 5 in Algorithm 4 is satisfied. Moreover, if $0\in\mathcal{S}$ and $y_{(x_{0},\theta_{0})}$ and its derivatives are known at $t=0$ , then we could choose $\tilde{t}=0$ if Step 5 in Algorithm 4 is satisfied and not need to integrate backwards.

Remark 5.

In the algorithms we present here, one can choose a suitable set of time instants $\{t_{j,1},\dots,t_{j,q_{j}}\}$ , $j\in\{1,\dots,p\}$ , for Algorithms 3 and 3, and almost any time $\tilde{t}$ for Algorithm 4 that fulfill the required properties to reconstruct univocally the parameters and initial condition. Then, one can integrate the system up to the current time $t$ to reconstruct the current state. Reconstructing parameters and initial condition or parameters and current state are equivalent as already mentioned, provided the measurements to be error free. However, when facing real data corrupted by some noise, the estimation might depend on the choice of the time instants $\{t_{j,\ell}\}_{\ell=1}^{q_{j}}$ , for $j\in\{1,\dots,p\}$ , or the time $\tilde{t}$ , and is tainted by some error. Then, integrating the system backwards up to time $0$ and forward up to time $t$ might propagate and amplify the estimation error on the initial condition and the current state, respectively. A way to improve this estimation is to use a filter in the spirit of an observer (see e.g. [37]). In this work, as the objective is to analyze the identifiability and observability properties, we do not consider here the robustness issue with respect to corrupted data. This will be the matter of a future work.

Up to now, we have provided some hypotheses that ensure a system to be observable, identifiable or jointly observable-identifiable, and some constructive algorithms to recover the initial condition and/or the parameter vector. Besides, if we have some analyticity properties, we can use almost any time $t\in\mathcal{I}$ to recover these unknowns. Epidemiological models based on autonomous ODEs are typically analytic systems. In particular, if $f$ is a combination of linear and bilinear terms of the state variables, then $f$ is clearly analytic in all $\mathbb{R}^{n}$ . Due to the Cauchy-Kovalevskaya theorem for ODEs (see [39]), each initial condition provides a unique locally analytic solution. In an autonomous system, any point of a solution can be considered as the initial condition of the same solution, and hence this solution is analytic everywhere. Moreover, in the example presented in Section 2, the output is $y=kI$ , which, since $I$ is analytic, is also analytic.

Remark 6.

We can apply the proposed algorithms to some systems with piece-wise constant parameters. If the parameter vector is $\theta_{s}$ in an interval $[t_{s},t_{s+1})$ , $s\in\{0,1,2,\dots\}$ , and we know $t_{s}$ , for all $s$ , then we can check the different hypotheses exposed in Algorithms 3, 3 and 4 considering for each $s$ that we know $y_{(x(t_{s}),\theta_{s})}$ in $[t_{s},t_{s+1})$ , and apply them respectively if the hypotheses are fulfilled.

Let us now illustrate all the presented theory through the classical SIRS model.

4 Application to the SIRS model

In this section, we study the observability, identifiability and joined observability-identifiability of the SIRS model (1), previously presented in Section 2, under the observation of an unknown fraction $k\in(0,1]$ of individuals. The associated ODE system is conservative, i.e., $\dot{S}+\dot{I}+\dot{R}=0$ ; in particular, $S$ is the fraction of susceptible individuals, $I$ is the fraction of infectious individuals and $R$ is the fraction of recovered individuals, and then $S+I+R=1$ , for all $t\geq 0$ (see, for example, [9, Chapter 10.4]). Then, we consider the following reduced system:

\left\{\begin{array}[]{ccl}\dot{S}&=&-\beta SI+\mu(1-S-I),\\ \dot{I}&=&\beta SI-\gamma I,\\ y&=&kI,\end{array}\quad\forall\,t\geq 0,\right.

(12)

for some initial condition $(S(0),I(0))^{\mathrm{T}}=(S_{0},I_{0})^{\mathrm{T}}\in\Omega=\{(\xi_{1},\xi_{2})^{\mathrm{T}}\in[0,1]^{2}\ :\ \xi_{1}+\xi_{2}\leq 1\}$ . We assume that the four parameters, $k$ , $\beta$ , $\gamma$ and $\mu$ , and/or the initial condition are unknown.

In this system, $\beta>0$ (days^-1) is the disease contact rate, $\gamma>0$ (days^-1) is the transition rate from compartment $I$ to compartment $R$ , and $\mu>0$ (days^-1) is the transition rate from compartment $R$ to compartment $S$ . The following results are well-known ([9, Chapter 10.4], [41]):

•

The system of ODEs of System (12) with initial condition in $\Omega$ has a unique solution.
•

The set $\Omega$ is positively invariant with respect to the system of ODEs of System (12).
•

The basic reproduction number is defined as $\mathcal{R}_{0}=\beta/\gamma$ . It is an approximation of the number of cases that one infected person generates on average over the course of its infectious period, in an uninfected population and without special control measures.
•

The (unique) disease-free equilibrium (DFE), which is $P_{\mathrm{f}}=(1,0)^{\mathrm{T}},$ is globally asymptotically stable when $\mathcal{R}_{0}<1$ and unstable when $\mathcal{R}_{0}>1$ ; in particular, in this case, it is a saddle point whose stable manifold is $I\equiv 0$ .
•

When $\mathcal{R}_{0}>1$ , there exists an endemic equilibrium (EE) in $\Omega$ given by $P_{\mathrm{e}}(\theta)=\left(\dfrac{1}{\mathcal{R}_{0}},\right.$ $\left.\mu\dfrac{1-1/\mathcal{R}_{0}}{\gamma+\mu}\right)^{\mathrm{T}},$ which is globally asymptotically stable out of $I\equiv 0$ .

Having done this quick wrap-up on the main characteristics of the SIRS model, we may start to study if it is observable, identifiable or jointly observable-identifiable. Notice that, in particular, the solutions of this system are analytic and, hence, so it is the output.

4.1 Identifiability and observability

Notice that the solutions and observations of System (12) are analytic in $\mathcal{I}$ . We will see that the procedure described in Algorithm 3 is appropriate for this case, and gives more general results than Algorithm 3. Hence, in the following, we are going to focus on Algorithm 3 considering we know $y_{(x_{0},\theta_{0})}$ in $[a,b)\subset\mathcal{I}$ . In particular, to check the observability, identifiability or joined observability-identifiability of System (12), we are going to follow Steps 1 and 2 and see if the required conditions can be satisfied. First of all, notice that the number of state variables of System (12) is $n=2$ , the number of unknown parameters is $b=4$ and the number of observations is $m=1$ . Let $\xi=(S_{0},I_{0})^{\mathrm{T}}\in\Omega$ , $\theta=(k,\beta,\gamma,\mu)^{\mathrm{T}}\in\Theta=(0,1]\times(0,\infty)^{3}$ , $x=(S,I)^{\mathrm{T}}$ . Then, System (12) can be consider as a system of the form given in (2), with

f(x,\theta)=\left(\begin{array}[]{c}-\beta x_{1}x_{2}+\mu(1-x_{1}-x_{2})\\ \beta x_{1}x_{2}-\gamma x_{2}\end{array}\right),\quad h(x,\theta)=kx_{2}.

In the following, we will find different suitable sets $\Omega_{1}$ and $\Omega_{2,\theta}$ , for any $\theta\in\Theta$ , for the above mentioned Steps 1 and 2, respectively. We will do it constructively.

Step 1. In this case, since $m=1$ , we have to find $d_{1}\in\mathbb{N}_{0}$ such that, for any $\theta\in\Theta$ , the following function is injective in some suitable set $\Omega_{1}\subset\Omega$ positively invariant w.r.t. the ODE system given in (12):

\mathcal{L}_{\theta,d_{1}}:\xi\mapsto\left(y_{(\xi,\theta)}^{(0)}(0),\dots,y_{(\xi,\theta)}^{(d_{1})}(0)\right).

For any $(\xi,\theta)\in\Omega\times\Theta$ , letting $x(t;\xi,\theta)=(S(t;\xi,\theta),I(t;\xi,\theta))^{\mathrm{T}}$ be the solution to System (12) with initial condition $\xi$ and parameter vector $\theta$ , we have

y^{(0)}_{(\xi,\theta)}(t)=h_{\theta}(x(t;\xi,\theta))=kI(t;\xi,\theta),\quad\forall\,t\geq 0.

Notice that $\mathcal{L}_{\theta,0}(\xi)=y^{(0)}_{(\xi,\theta)}(0)=k\xi_{2}$ is injective in a set $\Omega_{1}\subset\Omega$ if, and only if, there exist $\Xi\subset[0,1]$ and some function $G:\Xi\rightarrow[0,1]$ such that $\Omega_{1}=\{(G(\xi_{2}),\xi_{2})\in\Omega\ :\ \xi_{2}\in\Xi\}$ . A set $\Omega_{1}$ that has this form is positively invariant with respect to the system of ODEs of System (12) if, and only if, for any initial condition $\xi\in\Omega_{1}$ , $S(t;\xi,\theta)=G(I(t;\xi,\theta))$ , for all $t\geq 0$ . Two examples of such orbits are $(S,I)\equiv P_{\mathrm{f}}$ and, whenever $\mathcal{R}_{0}>1$ , $(S,I)\equiv P_{\mathrm{e}}(\theta)$ .

These type of sets are very restrictive and not very useful in applications. Therefore, let us consider $d_{1}=1$ instead of $d_{1}=0$ . We will show that this is a more convenient choice since the positively invariant set w.r.t. the ODE system given in (12) that we will derive from it is less restrictive.

Differentiating $y_{(\xi,\theta)}^{(0)}(t)$ one time leads to the following equation:

\begin{array}[]{ll}y^{(1)}_{(\xi,\theta)}(t)&=\mathrm{D}h_{\theta}(x(t;\xi,\theta))f_{\theta}(x(t;\xi,\theta))=(0,k)f_{\theta}(x(t;\xi,\theta))\\ &=k(\beta S(t;\xi,\theta)-\gamma)I(t;\xi,\theta)(t)=(\beta S(t;\xi,\theta)-\gamma)y^{(0)}_{(\xi,\theta)}(t),\quad\forall\,t\geq 0.\end{array}

Let us consider $\Omega_{1}=\{(\xi_{1},\xi_{2})^{\mathrm{T}}\in\Omega\ :\ \xi_{2}\neq 0\}=\{(\xi_{1},\xi_{2})^{\mathrm{T}}\in[0,1)\times(0,1]\ :\ \xi_{1}+\xi_{2}\leq 1\}$ .

Lemma 5.

For any $\theta=(k,\beta,\gamma,\mu)^{\mathrm{T}}\in\Theta$ , the following map is injective:

\mathcal{L}_{\theta,1}:\Omega_{1}\longrightarrow\mathbb{R}^{2}:\xi\longmapsto\left(y^{(0)}_{(\xi,\theta)}(0),y^{(1)}_{(\xi,\theta)}(0)\right).

(13)

Proof.

The proof is straightforward, since, given $\theta=(k,\beta,\gamma,\mu)\in\Theta$ , $\xi\in\Omega_{1}$ ,

y^{(0)}_{(\xi,\theta)}(0)=k\xi_{2},\quad y^{(1)}_{(\xi,\theta)}(0)=k(\beta\xi_{1}-\gamma)\xi_{2}.

Therefore, it is easy to see that (13) is injective. ∎

Notice that, if we change $\Omega_{1}$ by $\Omega$ in Lemma 5, then the mapping $\mathcal{L}_{\theta,1}$ is not injective.

To continue with Step 1 of Algorithm 3, we point out that $\Omega_{1}$ is positively invariant with respect to the system of ODEs of System (12). Notice that $\Omega_{1}$ is the same set as $\Omega$ except for the manifold $\{(\xi_{1},0):\xi_{1}\in[0,1]\}$ , which is also positively invariant w.r.t. the ODE system given in (12) (recall that this is the stable manifold associated to $P_{\mathrm{f}}$ when $\mathcal{R}_{0}>1$ ). Then, due to the uniqueness of solutions, $\Omega_{1}$ is also positively invariant w.r.t. the ODE system given in (12). Hence, System (12) is observable on $\Omega_{1}$ in any semi-open interval $[a,b)\subset\mathcal{I}$ with parameters in $\Theta$ . If we know $\theta_{0}\in\Theta$ and some observations $y_{(x_{0},\theta_{0})}(t)$ , for all $t\in[a,b)$ , then we can determine $x_{0}$ performing now Steps 3, 7 and 8 in Algorithm 3.

Step 2. We now need to find a suitable sets $\Omega_{2,\theta}\subset\Omega$ , for every $\theta\in\Theta$ , positively invariant w.r.t. the ODE system given in (12), and maps $g:\mathcal{D}\rightarrow\mathbb{R}^{q+p}$ and $r:\Theta\rightarrow\mathbb{R}^{p}$ , for some suitable $\mathcal{D}\subset\mathbb{R}^{d_{1}^{\prime}+1}$ , $d_{1}^{\prime}\in\mathbb{N}\cup\{0\}$ , $q,p\in\mathbb{N}$ , such that (C1), (C2) and (C3) of Theorem 2 are satisfied for some connected $\mathcal{S}\subset\mathcal{I}$ such that $[a,b)\subset\mathcal{S}$ and $g_{j,l}\left(y_{(\xi,\theta)}^{(0)}(t),\dots,y_{(\xi,\theta)}^{(d_{1}^{\prime})}(t)\right)$ are analytic w.r.t. $t\in\mathcal{S}$ , for any $j\in\{1,\dots,p\}$ , $l\in\{1,\dots,q_{j}\}$ , and $(\xi,\theta)\in\Gamma_{2,\Theta}$ .

In the following, for the sake of cleanliness, we make a slight abuse of notation and avoid the specification of $\xi$ , $\theta$ and $t$ when no confusion is possible; besides, we denote $y^{(0)}$ and $y^{(1)}$ as $y$ and $\dot{y}$ , respectively.

Let us start by checking (C1) of Theorem 2. Since

y=kI,\quad\dot{y}=(\beta S-\gamma)y,\quad\forall\ t\geq 0,

(14)

neither $d_{1}^{\prime}=0$ nor $d_{1}^{\prime}=1$ are suitable. In particular, we cannot obtain a function $r$ injective in $\Theta$ suitable for (7) since we do not have any information on $\mu$ in the expressions of $y^{(0)}$ and $y^{(1)}$ given in (14).

Then, if we continue differentiating $\dot{y}$ , denoting $y^{(2)}$ as $\ddot{y}$ , we obtain:

\ddot{y}=(\beta S-\gamma)\dot{y}+\beta\dot{S}y=(\beta S-\gamma)y+\beta(-\beta SI+\mu(1-S-I))y,\quad\forall\ t\geq 0.

(15)

Notice that, if $y\not\equiv 0$ , which holds when we consider initial conditions in $\Omega_{1}$ , we can define $(S,I)$ in terms of $y$ , $\dot{y}$ and $\theta$ from (14) as

S=\dfrac{1}{\beta}\left(\dfrac{\dot{y}}{y}+\gamma\right),\quad I=\dfrac{y}{k},\quad\forall\ t\geq 0.

(16)

Hence, substituting (16) in (15) and performing some computations, we obtain an equation in the form of (7) in (C1) of Theorem 2:

\ddot{y}-\dfrac{\dot{y}^{2}}{y}=-y\mu(\gamma-\beta)-y^{2}\dfrac{\beta}{k}(\gamma+\mu)-\dot{y}\mu-y\dot{y}\dfrac{\beta}{k},\quad\forall\ t\geq 0.

(17)

Then, we can see that (C1) of Theorem 2 is satisfied with $d_{1}^{\prime}=2$ , $p=1$ , $q_{1}=q=4$ ,

g_{0}(y,\dot{y},\ddot{y})=\ddot{y}-\dfrac{\dot{y}^{2}}{y},\ g_{1}(y,\dot{y},\ddot{y})=-y,\ g_{2}(y,\dot{y},\ddot{y})=-y^{2},\ g_{3}(y,\dot{y},\ddot{y})=-\dot{y},\ g_{4}(y,\dot{y},\ddot{y})=-y\dot{y},

and

r(\theta)=\left(\mu(\gamma-\beta),\dfrac{\beta}{k}(\gamma+\mu),\mu,\dfrac{\beta}{k}\right),\quad\forall\,\theta\in\Theta.

Since $p=1$ , we have denoted $g_{1,l}=g_{l}$ , $l\in\{0,\dots,4\}$ , for the sake of simplicity.

Let us now check (C2) of Theorem 2, i.e., if the function $r:\Theta\rightarrow r(\Theta)\subset\mathbb{R}^{4}$ is injective. Indeed, given some parameter vector $\theta=(k,\beta,\gamma,\mu)^{\mathrm{T}}\in\Theta$ , it is easy to see that we can invert the equation $r(\theta)=\sigma$ and obtain univocally

\left(k,\beta,\gamma,\mu\right)=\left(\dfrac{1}{\sigma_{4}}\left(\dfrac{\sigma_{2}}{\sigma_{4}}-\sigma_{3}-\dfrac{\sigma_{1}}{\sigma_{3}}\right),\dfrac{\sigma_{2}}{\sigma_{4}}-\sigma_{3}-\dfrac{\sigma_{1}}{\sigma_{3}},\dfrac{\sigma_{2}}{\sigma_{4}}-\sigma_{3},\sigma_{3}\right).

Finally, we need to check (C3) of Theorem 2, i.e., if, for any $(\xi,\theta)\in\Gamma_{2,\Theta}$ , the functions $\{g_{i}(y,\dot{y},\ddot{y})\}_{i=1}^{4}$ are linearly independent with respect to time in some suitable $\mathcal{S}\subset\mathcal{I}$ , for initial conditions in some suitable sets $\Omega_{2,\theta}\subset\Omega$ , for $\theta\in\Theta$ , positively invariant w.r.t. the ODE system given in (12). Since these functions are defined in all $\mathcal{I}$ , we will consider $\mathcal{S}=\mathcal{I}$ . Then, if they are linearly independent in $\mathcal{I}$ , given that they are clearly also analytic, we will have linear independence in any $[a,b)\subset\mathcal{I}$ . In the following Lemma 6, we will directly prove linear independence for any $[a,b)\subset\mathcal{I}$ .

Lemma 6.

For any $[a,b)\subset\mathcal{I}$ , the functions $\{g_{i}(y,\dot{y},\ddot{y})\}_{i=1}^{4}$ are linearly independent in $[a,b)$ if $\Gamma_{2,\Theta}=\bigcup_{\theta\in\Theta}(\Omega_{1}\setminus\{P_{\mathrm{e}}(\theta)\})\times\{\theta\}$ and $(\xi,\theta)\in\Gamma_{2,\Theta}$ , i.e., $\Omega_{2,\theta}=\Omega_{1}\setminus\{P_{\mathrm{e}}(\theta)\}$ , where $P_{\mathrm{e}}(\theta)$ denotes the endemic equilibrium point associated with $\theta$ .

Proof.

Consider some $[a,b)\subset\mathcal{I}$ . For some $\theta\in\Theta$ , we consider $\xi\in\Omega$ , and assume there exist $a_{1},a_{2},a_{3},a_{4}\in\mathbb{R}$ such that

a_{1}y(t)+a_{2}y(t)^{2}+a_{3}\dot{y}(t)+a_{4}y(t)\dot{y}(t)=0,\quad\forall\ t\in[a,b).

(18)

We will prove that $a_{1}=a_{2}=a_{3}=a_{4}=0$ if $\xi\in\Omega_{2,\theta}$ , for some $\Omega_{2,\theta}\subset\Omega$ . As mentioned in Step 1, if $\xi_{2}=0$ , then $y\equiv 0$ in $\mathcal{I}$ , and hence (18) is true for values $a_{1},a_{2},a_{3},a_{4}$ not necessarily satisfying $a_{1}=a_{2}=a_{3}=a_{4}=0$ , which implies that the corresponding functions $\{g_{i}(y,\dot{y},\ddot{y})\}_{i=1}^{4}$ are not linearly independent in $[a,b)$ .

Therefore, let us consider $\xi_{2}\neq 0$ , i.e., $\xi\in\Omega_{1}$ . In this case, if $(S(t),I(t))^{\mathrm{T}}$ is the solution of the ODE system of (12) with initial condition $\xi$ and parameter vector $\theta$ , recall that $\Omega_{1}$ is positively invariant w.r.t. the ODE system given in (12), i.e., $(S(t),I(t))^{\mathrm{T}}\in\Omega_{1}$ , for all $t\geq 0$ . Then, we can rewrite (18) as

	$\displaystyle a_{1}kI(t)+a_{2}k^{2}I^{2}(t)+a_{3}kI(t)(\beta S(t)-\gamma)+a_{4}k^{2}I^{2}(t)(\beta S(t)-\gamma)$	$\displaystyle=$	$\displaystyle 0\implies$
	$\displaystyle\overset{(S,I)^{\mathrm{T}}\in\Omega_{1}}{\implies}a_{1}+a_{2}kI(t)+a_{3}(\beta S(t)-\gamma)+a_{4}kI(t)(\beta S(t)-\gamma)$	$\displaystyle=$	$\displaystyle 0,\quad\forall\,t\in[a,b),$

which leads to

(a_{1}-a_{3}\gamma)+a_{3}\beta S(t)+(a_{2}k-a_{4}k\gamma)I(t)+a_{4}k\beta S(t)I(t)=0,\quad\forall\,t\in[a,b).

Determining $\{a_{i}\}_{i=1}^{4}$ not all of them null such that this is fulfilled is equivalent to determining $A,B,C,D$ not all of them null such that

A+BS(t)+CI(t)+DS(t)I(t)=0,\quad t\in[a,b),

(19)

since $A=B=C=D=0$ if, and only if, $a_{1}=a_{2}=a_{3}=a_{4}=0$ . That is, we will check whether $1,S,I,SI$ are linearly independent functions in $[a,b)$ considering initial conditions in $\Omega_{1}$ or not.

If we consider that $S$ or $I$ are constant, (19) is true for values of $A,B,C,D$ not necessarily satisfying $A=B=C=D=0$ , which implies that the corresponding functions $\{g_{i}(y,\dot{y},\ddot{y})\}_{i=1}^{4}$ are not linearly independent in $[a,b)$ . In $\Omega_{1}$ , due to the analyticity of the system, this is only possible if $(S,I)^{\mathrm{T}}\equiv P_{\mathrm{e}}(\theta)$ . Indeed, $\dot{I}\equiv 0$ implies $I\equiv 0$ , which is excluded from $\Omega_{1}$ , or $S\equiv\beta/\gamma$ and hence $(S,I)^{\mathrm{T}}\equiv P_{\mathrm{e}}(\theta)$ ; on the other hand, if $\dot{S}\equiv 0$ , then, $S\equiv\xi_{1}$ and

I\equiv\dfrac{\mu(1-\xi_{1})}{\mu+\beta\xi_{1}},\quad\forall\,t\in[a,b),

i.e., we are at an equilibrium point.

Hence, let us take $\Omega_{2,\theta}=\Omega_{1}\setminus\{P_{\mathrm{e}}(\theta)\}$ and check if, when $\xi\in\Omega_{2,\theta}$ , the corresponding functions $1,S,I,SI$ are linearly independent in $\Omega_{2,\theta}$ for times in $[a,b)$ . Assume that (19) holds and at least one between $C$ and $D$ is non-null. If we consider initial conditions in $\Omega_{2,\theta}$ , we obtain the following expression for $I$ :

I(t)=-\dfrac{A+BS(t)}{C+DS(t)},

where $C+DS(t)$ is non-null for almost every $t\in[a,b)$ , since $S$ is analytic and non-constant in $\Omega_{2,\theta}$ . Without loss of generality, assume $C+DS(t)$ is non-null for all $t\in[a,b)$ . We can hence differentiate this expression for $I(t)$ and obtain

\dot{I}(t)=-\dot{S}(t)\dfrac{BC-AD}{(C+DS(t))^{2}}=\left((\beta S(t)+\mu)\dfrac{A+BS(t)}{C+DS(t)}-\mu(1-S(t))\right)\dfrac{BC-AD}{(C+DS(t))^{2}},\quad\forall\,t\in[a,b).

On the other hand,

\dot{I}(t)=\beta S(t)I(t)-\gamma I(t)=(\gamma-\beta S(t))\dfrac{A+BS(t)}{C+DS(t)},\quad\forall\,t\in[a,b).

If we equal both expressions for $\dot{I}(t)$ , we get

\left((\beta S(t)+\mu)\dfrac{A+BS(t)}{C+DS(t)}-\mu(1-S(t))\right)\dfrac{BC-AD}{(C+DS(t))^{2}}=(\gamma-\beta S(t))\dfrac{A+BS(t)}{C+DS(t)},

for $t\in[a,b)$ . After some computations in order to get rid of the denominators, we reach the following polynomial on $S(t)$ :

c_{4}S(t)^{4}+c_{3}S(t)^{3}+c_{2}S(t)^{2}+c_{1}S(t)+c_{0}=0,\quad\forall\,t\in[a,b),

where

$\displaystyle c_{0}$	$\displaystyle=$	$\displaystyle\mu ABC-\mu A^{2}D+\mu ACD-\gamma AC^{2}-\mu BC^{2},$
$\displaystyle c_{1}$	$\displaystyle=$	$\displaystyle\mu AD^{2}-\beta A^{2}D-(2\gamma+\mu)ACD-\mu ABD-\mu BCD+\beta AC^{2}-(\gamma-\mu)BC^{2}+\mu B^{2}C+\beta ABC,$
$\displaystyle c_{2}$	$\displaystyle=$	$\displaystyle\beta BC^{2}-(2\gamma-\mu)BCD-(\gamma+\mu)AD^{2}-\beta ABD+2\beta ACD+\beta B^{2}C,$
$\displaystyle c_{3}$	$\displaystyle=$	$\displaystyle\beta AD^{2}-\gamma BD^{2}+2\beta BCD,$
$\displaystyle c_{4}$	$\displaystyle=$	$\displaystyle\beta BD^{2}.$

Since $S$ is analytic non-constant in $[a,b)$ , for this polynomial to be 0 in $[a,b)$ , we need that $c_{0}=\dots=c_{4}=0$ . We have assumed that at least one between $C$ and $D$ is non-null. We have two subcases:

•

Assume $D\neq 0$ . Then, it is straightforward observing that $c_{3}=c_{4}=0$ implies $A=B=0$ . But then, this implies that $I\equiv 0$ in $[a,b)$ , which cannot happen for solutions in $\Omega_{2,\theta}$ . Hence, $D=0$ .
•

Assume $C\neq 0$ , and $D=0$ due to the previous argument. Then, $c_{2}=0$ implies $B=0$ or $B=-C$ . If $B=0$ , from $c_{1}=0$ , we obtain $A=0$ , which again cannot happen for solutions in $\Omega_{2,\theta}$ . If we consider $B=-C$ , $c_{1}=0$ implies $\gamma C^{3}=0$ , which is not true since $\gamma>0$ and we have assumed $C\neq 0$ .

Hence, we must have $C=D=0$ . Then, equation (19) is

A+BS(t)=0,\quad\forall\,t\in[a,b).

But $S$ is analytic non-constant for solutions in $\Omega_{2,\theta}$ , and hence $A=B=0$ .

Therefore, for any $\theta\in\Theta$ , $1,S,I,SI$ or, equivalently, $\{g_{i}(y,\dot{y},\ddot{y})\}_{i=1}^{4}$ are linearly independent in any semi-open subinterval $[a,b)\subset\mathcal{I}$ if we consider $\xi$ in the following set:

\Omega_{2,\theta}=\left\{(\xi_{1},\xi_{2})^{\mathrm{T}}\in[0,1)\times(0,1]\ :\ \xi_{1}+\xi_{2}\leq 1\right\}\setminus\{P_{\mathrm{e}}(\theta)\}\subset\Omega_{1}\subset\Omega.

∎

To continue with Step 2, we need to check if the set $\Omega_{2,\theta}$ is positively invariant with respect to the system of ODEs of System (12).

Lemma 7.

Let $\theta\in\Theta$ . The set $\Omega_{2,\theta}=\{(\xi_{1},\xi_{2})^{\mathrm{T}}\in[0,1)\times(0,1]\ :\ \xi_{1}+\xi_{2}\leq 1\}\setminus\{P_{\mathrm{e}}(\theta)\}$ is positively invariant with respect to the system

\left\{\begin{array}[]{ccl}\dot{S}&=&-\beta SI+\mu(1-S-I),\\ \dot{I}&=&\beta SI-\gamma I.\end{array}\right.

(20)

Proof.

As said at the beginning of Section 4 the set

\Omega=\{(\xi_{1},\xi_{2})^{\mathrm{T}}\in[0,1]^{2}\,:\,\xi_{1}+\xi_{2}\leq 1\}

is positively invariant with respect to System (20).

Let $P_{\mathrm{f}}=(1,0)^{\mathrm{T}}$ be the disease-free equilibrium of System (20) and, given $\theta\in\Theta$ , $P_{\mathrm{e}}(\theta)=(1/\mathcal{R}_{0},\mu(1-1/\mathcal{R}_{0})/(\mu+\gamma))^{\mathrm{T}}$ the endemic equilibrium (which is inside $\Omega$ if, and only if, $\mathcal{R}_{0}>1$ ). Since the solution of this system given an initial condition is unique in $t\geq 0$ , then the set $\Omega\setminus\{P_{\mathrm{f}},P_{\mathrm{e}}(\theta)\}$ is still positively invariant with respect to System (20).

On the other hand, if an initial condition $\xi$ is such that $\xi_{2}\neq 0$ , then $I(t)>0$ , for all $t\geq 0$ , due to the uniqueness of solutions. Indeed, since $S(t)\geq 0$ , $\forall\,t\geq 0$ ,

\dot{I}(t)=(\beta S(t)-\gamma)I(t)\geq-\gamma I(t),\quad\forall\,t\geq 0\quad\implies\quad I(t)\geq\xi_{2}\mathrm{e}^{-\gamma t}>0,\quad\forall\,t\geq 0.

Thus, $\{(\xi_{1},\xi_{2})^{\mathrm{T}}\in[0,1)\times(0,1]\ :\ \xi_{1}+\xi_{2}\leq 1\}\setminus\{P_{\mathrm{e}}(\theta)\}$ is also a positively invariant set with respect to System (20). ∎

Hence, Step 2 is finished, which allows us to conclude that System (12) is identifiable on $\Theta$ in any semi-open interval $[a,b)\subset\mathcal{I}$ with initial conditions consistent with the family $\{\Omega_{2,\theta}\}_{\theta\in\Theta}$ . Actually, if we know $x_{0}\in\Omega_{2,\theta_{0}}$ and some observations $y_{(x_{0},\theta_{0})}(t)$ , for $t\in[a,b)\subset\mathcal{I}$ , then we can determine $\theta_{0}$ performing Steps 3–6 in Algorithm 3.

Step 3. Set $\Omega_{\theta}=\Omega_{1}\cap\Omega_{2,\theta}=\Omega_{2,\theta}$ , $\theta\in\Theta$ . Then, this is a positively invariant set, and we assume $x_{0}\in\Omega_{\theta_{0}}$ .

Notice that, for any $\theta\in\Theta$ , not considering initial points in $\{(\xi_{1},0):\xi_{1}\in[0,1]\}\cup\{P_{\mathrm{e}}(\theta)\}$ is not restrictive with respect to the observations we can work with, since $y\equiv 0$ would mean there is no disease and therefore no point of study, and $y\equiv$ constant $\neq 0$ means the disease is already endemic and it will not vary unless we perturb the system (e.g., with migration or vaccination).

Hence, System (12) is jointly observable-identifiable on $\Gamma_{\Theta}$ in any $[a,b)\subset\mathcal{I}$ . If we had some particular observations $y_{(x_{0},\theta_{0})}$ in $[a,b)\subset\mathcal{I}$ , we could continue with Steps 4–8 in Algorithm 3 or Steps 4–9 in Algorithm 4.

4.2 A limiting case: The SIR model

Another basic compartmental model is the SIR model, where the considered populations are the same as for System (12), but there is no loss of immunity, i.e., $\mu=0$ . Again, we can consider the observation of an unknown portion of infected individuals and hence obtain the following system:

\left\{\begin{array}[]{lcl}\dot{S}&=&-\beta SI,\\ \dot{I}&=&\beta SI-\gamma I,\\ y&=&kI,\end{array}\right.

(21)

where we have already simplified the dimension taking into account that the quantity $S+I+R=1$ is preserved.

If one thinks of whether System (21) is observable, identifiable or jointly observable-identifiable, it seems it will satisfy the same properties as the previous System (12), and it will be easier to prove them, since we got rid of parameter $\mu$ . However, this is far from being true.

If we perform the same computations as for the SIRS model, we obtain the same result for Step 1, which implies that System (21) is observable on $\Omega_{1}=\{(\xi_{1},\xi_{2})^{\mathrm{T}}\in[0,1)\times(0,1]\ :\ \xi_{1}+\xi_{2}\leq 1\}$ in any $[a,b)\subset\mathcal{I}$ with parameters in $\Theta=(0,1]\times(0,\infty)^{2}$ . Regarding Step 2, as in (16) and (17) with $\mu=0$ , we obtain that

I=\dfrac{y}{k},\quad S=\dfrac{1}{\beta}\left(\dfrac{\dot{y}}{y}+\gamma\right)

(22)

and

\ddot{y}-\dfrac{\dot{y}}{y}=-\dfrac{\beta\gamma}{k}y^{2}-\dfrac{\beta}{k}y\dot{y}.

(23)

Then,

r(k,\beta,\gamma)=\left(\dfrac{\beta\gamma}{k},\dfrac{\beta}{k}\right),

which is clearly not injective in $\Theta$ , and hence we cannot complete Step 2 using (23). However, recalling Remark 1, System (21) is indeed identifiable considering a suitable positively invariant subset of $\Omega_{1}$ , as we can see in the following Lemma (8), which does not need the injectivity of $r$ .

Lemma 8.

System (21) is identifiable on $\Theta=(0,1]\times(0,\infty)^{2}$ in any $[a,b)\subset\mathcal{I}$ , with initial conditions in $\tilde{\Omega}=\{(\xi_{1},\xi_{2})^{\mathrm{T}}\in(0,1)^{2}\ :\ \xi_{1}+\xi_{2}\leq 1\}$ ²²2Note that $\tilde{\Omega}$ is independent of the choice of $\theta$ ..

Proof.

System (21) will be identifiable on $\Theta$ with initial conditions in $\tilde{\Omega}$ if, given any initial condition $\xi\in\tilde{\Omega}$ and some observations $y$ in $[a,b)$ , we can determine $\theta\in\Theta$ univocally.

First of all, one can prove, similarly to Lemma 7, that $\tilde{\Omega}$ is positively invariant with respect to the ODE system in System (21). Equations (22) and (23) are still valid. Therefore, if we prove that $y^{2}$ and $y\dot{y}$ are linearly independent w.r.t. $t\in[a,b)$ , then there is a unique solution $\sigma=(\sigma_{1},\sigma_{2})$ to

\ddot{y}-\dfrac{\dot{y}}{y}=-\sigma_{1}y^{2}-\sigma_{2}y\dot{y},

(24)

where $\sigma_{1}=\beta\gamma/k$ and $\sigma_{2}=\beta/k$ . To see that the linear independence holds, let us consider $a_{1},a_{2}\in\mathbb{R}$ such that

a_{1}y^{2}+a_{2}y\dot{y}=0,\quad\forall\ t\in[a,b).

Since $\tilde{\Omega}$ is positively invariant, then $I\not\equiv 0$ in $[a,b)$ , and hence $y\not\equiv 0$ , which implies

a_{1}y+a_{2}\dot{y}=0.

This is equivalent to

a_{1}kI+a_{2}kI(\beta S-\gamma)=0,\ \forall\ t\in[a,b),\quad\text{which implies}\quad a_{1}+a_{2}(\beta S-\gamma)=0,\ \forall\ t\in[a,b).

Since $S$ is analytic, $a_{1}$ or $a_{2}$ are non-null if, and only if, $S$ is constant. Given that $\dot{S}=-\beta SI$ , this can only happen if $S(t)=0$ or $I(t)=0$ , for all $t\geq 0$ , which cannot occur in $\tilde{\Omega}$ . Hence, $y^{2}$ and $y\dot{y}$ are linearly independent when we consider $\xi\in\tilde{\Omega}$ and $\theta\in\Theta$ , and therefore there is a unique solution $\sigma$ to (24). Then, we can determine univocally $\beta/k=\sigma_{2}$ and $\gamma=\sigma_{1}/\sigma_{2}$ .

If $a=0$ , taking into account that $\xi_{1},\xi_{2}\neq 0$ in $\tilde{\Omega}$ , we can conclude directly determining

k=\dfrac{y(0)}{\xi_{2}},\quad\gamma=\dfrac{\sigma_{1}}{\sigma_{2}},\quad\beta=\dfrac{1}{\xi_{1}}\left(\dfrac{\dot{y}(0)}{y(0)}+\dfrac{\sigma_{1}}{\sigma_{2}}\right).

If $a>0$ , notice now that we can make the following change of variables in the ODE system in System (21): $X=\beta S$ and $Y=kI$ . Then,

\left\{\begin{array}[]{lcl}\dot{X}&=&-\sigma_{2}XY,\\ \dot{Y}&=&Y(X-\gamma),\end{array}\right.

and we can recover $\beta\xi_{1}$ and $k\xi_{2}$ integrating backwards considering as initial condition

X(a)=\dfrac{\dot{y}(a)}{y(a)}+\dfrac{\sigma_{1}}{\sigma_{2}},\quad Y(a)=y(a).

Once we know $X(0)=\beta\xi_{1}$ and $Y(0)=k\xi_{2}$ , taking into account that $\xi_{1},\xi_{2}\neq 0$ in $\tilde{\Omega}$ , we can determine univocally $k$ , $\beta$ and $\gamma$ as

k=\dfrac{Y(0)}{\xi_{2}},\quad\beta=\dfrac{X(0)}{\xi_{1}},\quad\gamma=\dfrac{\sigma_{1}}{\sigma_{2}}.

∎

Therefore, System (21) is observable on $\Omega_{1}$ with parameters in $\Theta$ and identifiable on $\Theta$ with initial conditions in $\tilde{\Omega}\subset\Omega_{1}$ , in any $[a,b)\subset\mathcal{I}$ . Nevertheless, it is not jointly observable-identifiable on $\tilde{\Omega}\times\Theta$ . Actually, in [13], the authors treat this case and prove that, assuming both $x_{0}=(S_{0},I_{0})\in\tilde{\Omega}$ and $\theta_{0}=(k_{0},\beta_{0},\gamma_{0})\in\Theta$ unknown, we can only determine $\gamma_{0}$ , $\beta_{0}/k_{0}$ , $\beta_{0}S_{0}$ and $k_{0}I_{0}$ , i.e., it is partially jointly observable-identifiable. In fact, considering (22) and (23), if $(\sigma_{1},\sigma_{2})$ is the solution to

\ddot{y}-\dfrac{\dot{y}}{y}=-\sigma_{1}y^{2}-\sigma_{2}y\dot{y}

(which we know is unique due to the proof of Lemma 8) and $(k_{0},\beta_{0},\gamma_{0})$ is any solution to $r(k_{0},\beta_{0},\gamma_{0})=(\sigma_{1},\sigma_{2})$ , notice that our method matches this result, since we have

I=\dfrac{y}{k_{0}},\quad S=\dfrac{1}{\beta_{0}}\left(\dfrac{\dot{y}}{y}+\gamma_{0}\right),\quad\dfrac{\beta_{0}\gamma_{0}}{k_{0}}=\sigma_{1},\quad\dfrac{\beta_{0}}{k_{0}}=\sigma_{2},

and hence, if we know $y(0)$ and $\dot{y}(0)$ , it is straightforward checking that we can recover $\gamma_{0}$ , $\beta_{0}/k_{0}$ , $\beta_{0}S_{0}$ , $k_{0}I_{0}$ univocally as follows:

\gamma_{0}=\dfrac{\sigma_{1}}{\sigma_{2}},\quad\dfrac{\beta_{0}}{k_{0}}=\sigma_{2},\quad\beta_{0}S_{0}=\dfrac{\dot{y}(0)}{y(0)}+\dfrac{\sigma_{1}}{\sigma_{2}},\quad k_{0}I_{0}=y(0),

and we cannot obtain more information. If $a>0$ , then we can proceed as in the proof of Lemma 8 performing the change of variables $X=\beta_{0}S$ and $Y=k_{0}I$ and integrating backwards from $t_{0}=a$ to $0$ .

One could also not know at first if some given observations $y_{(x_{0},\theta_{0})}$ correspond to an SIR or an SIRS model when both $x_{0}$ and $\theta_{0}$ are unknown, and wonder if the same data might be reproduced with two different sets of parameters $(k_{1},\beta_{1},\gamma_{1},0)^{\mathrm{T}}$ and $(k_{2},\beta_{2},\gamma_{2},\mu_{2})^{\mathrm{T}}$ , $\mu_{2}\neq 0$ . To tackle this question, we consider an extended SIRS model with parameters

\theta:=(k,\beta,\gamma,\mu)^{\mathrm{T}}\in\Theta^{*}=(0,1]\times(0,\infty)^{2}\times[0,\infty).

Notice that the SIR model is a particular case of this extension of the SIRS model, and can be regarded as a limiting case of the SIRS model when $\mu\rightarrow 0$ .

Then again, in Step 2 we would obtain

r(k,\beta,\gamma,\mu)=\left(\mu(\gamma-\beta),\dfrac{\beta}{k}(\gamma+\mu),\mu,\dfrac{\beta}{k}\right),

which is not injective in $\Theta^{*}$ and hence we cannot complete this step with this function $r$ . Nevertheless, this extended model can help to distinguish whether the observations come from an SIR model or an SIRS model, as illustrated in Section 4.2.2. Let us first do a quick comparison between both models in Section 4.2.1.

4.2.1 Comparison between SIR and SIRS models

Although passing from $\mu=0$ to $\mu>0$ may change substantially the behavior of the solutions, since the SIR model does not admit an endemic state, whereas the SIRS model does, they can be hardly distinguishable at early stages if $\mu$ is very small. We let $x_{\mathrm{SIR}}(t)$ and $x_{\mathrm{SIRS}}(t)$ be the solutions to the SIR given by the ODE system of (21) and the SIRS model (20), respectively, with the same initial condition, and we study the dependence of $||x_{\mathrm{SIR}}(t)-x_{\mathrm{SIRS}}(t)||_{2}$ on $\mu$ . To do this, we are going to consider the same parameters $\beta$ and $\gamma$ for both models, and we will base ourselves on Theorem 3.4., Chapter 3 of [41], particularized to our autonomous context:

Theorem 4.

[41, Chapter 3] Let $f$ be a Lipschitz map on $W$ with a global Lipschitz constant $L$ , where $W\subset\mathbb{R}^{n}$ is an open connected set. Let $y(t)$ and $z(t)$ be solutions of

\dot{y}=f(y),\ y(t_{0})=y_{0},\quad\text{and}\quad\dot{z}=f(z)+g(z),\ z(t_{0})=z_{0},

such that $y(t)$ , $z(t)\in W$ , for all $t\in[t_{0},t_{1}]$ . Suppose that

||g(x)||\leq\mu,\quad\forall\,x\in W,

for some $\mu>0$ . Then,

||y(t)-z(t)||\leq||y_{0}-z_{0}||\exp(L(t-t_{0}))+\dfrac{\mu}{L}[\exp(L(t-t_{0}))-1],\quad\forall\,t\in[t_{0},t_{1}].

Then, we are going to check the different conditions required in Theorem 4 in order to obtain an estimation of $||x_{\mathrm{SIR}}(t)-x_{\mathrm{SIRS}}(t)||_{2}$ .

Notice that both models share the same positively invariant set

\Omega=\{(\xi_{1},\xi_{2})^{\mathrm{T}}\in[0,1]^{2}\ :\ \xi_{1}+\xi_{2}\leq 1\},

which is compact, but the derivatives are well defined. Let $x_{\mathrm{SIR}}=(S_{\mathrm{SIR}},I_{\mathrm{SIR}})^{\mathrm{T}}$ and $x_{\mathrm{SIRS}}=(S_{\mathrm{SIRS}},I_{\mathrm{SIRS}})^{\mathrm{T}}$ be the state variables associated to the ODE system of (21) and (20), respectively. We are going to rewrite these systems in the form presented in Theorem 4. Let $x_{0}$ be the same initial condition at $t_{0}$ for both models, $\eta=(\beta,\gamma,\mu)^{\mathrm{T}}$ and $\nu=(\beta,\gamma)^{\mathrm{T}}$ two parameter vectors, and the following function:

f_{\mathrm{SIR}}(x,\nu)=\left(\begin{array}[]{c}-\beta x_{1}x_{2}\\ \beta x_{1}x_{2}-\gamma x_{2}\end{array}\right),\quad x\in\Omega.

Notice that $f_{\mathrm{SIR}}(x,\nu)$ is Lipschitz in $x$ on the compact set $\Omega$ , for some Lipschitz constant $L>0$ which is independent of $\mu$ .

Let us also define

g(x,\mu)=\left(\begin{array}[]{c}\mu(1-x_{1}-x_{2})\\ 0\end{array}\right),\quad x\in\Omega.

Then, $g$ satisfies

||g(x,\mu)||_{2}=\left|\left|\left(\begin{array}[]{c}\mu(1-x_{1}-x_{2})\\ 0\end{array}\right)\right|\right|_{2}\leq\mu,\quad\forall\,x\in\Omega.

We consider now the two following systems:

\left(\begin{array}[]{c}\dot{S}_{\mathrm{SIR}}(t;x_{0},\nu)\\ \dot{I}_{\mathrm{SIR}}(t;x_{0},\nu)\end{array}\right)=f_{\mathrm{SIR}}(x_{\mathrm{SIR}}(t;x_{0},\nu),\nu)

and

\left(\begin{array}[]{c}\dot{S}_{\mathrm{SIRS}}(t;x_{0},\eta)\\ \dot{I}_{\mathrm{SIRS}}(t;x_{0},\eta)\end{array}\right)=f_{\mathrm{SIR}}(x_{\mathrm{SIRS}}(t;x_{0},\eta),\nu)+g(x_{\mathrm{SIRS}}(t;x_{0},\eta),\mu).

Then, we are under the conditions of Theorem 4. Let us make an abuse of notation and let $x_{\mathrm{SIR}}(t)=x_{\mathrm{SIR}}(t;x_{0},\nu)$ and $x_{\mathrm{SIRS}}(t)=x_{\mathrm{SIRS}}(t;x_{0},\eta)$ . It is fulfilled that

||x_{\mathrm{SIR}}(t)-x_{\mathrm{SIRS}}(t)||_{2}\leq\dfrac{\mu}{L}[\exp(L(t-t_{0}))-1],\quad\forall\,t\geq t_{0}.

This is, for any $\varepsilon>0$ , $\bar{t}>t_{0}$ , there exists some $\mu>0$ small enough such that

||x_{\mathrm{SIR}}(t)-x_{\mathrm{SIRS}}(t)||_{2}\leq\dfrac{\mu}{L}[\exp(L(t-t_{0}))-1]<\varepsilon,\quad\forall\,t\in[t_{0},\bar{t}].

(25)

To illustrate this, we can observe in Figure 1 (Left) how both solutions are hardly distinguishable when considering $\beta=2.5$ , $\gamma=1$ , $\mu=0.001$ and $x_{\mathrm{SIR}}(0)=x_{\mathrm{SIRS}}(0)=(0.9,0.1)^{\mathrm{T}}$ . Moreover, the infectious compartment $I_{\mathrm{SIRS}}$ presents a slow-fast behavior for the SIRS model when near to the invariant manifold $I_{\mathrm{SIRS}}\equiv 0$ ; we observe in Figure 1 (Right, Bottom) how it takes a lot of time to move away from the manifold and then approaches it very fast, and hence remains most of the time very close to the solution of $I_{SIR}$ for the SIR model. Hence, it is reasonable that $\mu=0$ is hardly distinguishable from small values of $\mu$ knowing only the infectious compartment. Moreover, Figure 1 illustrates that this difficulty may not only be at initial times, but also for intermediate intervals of time.

Refer to caption — Figure 1: Comparison between the solutions of an SIR and an SIRS models with same initial condition $(0.9,0.1)$ and parameters $\beta=2.5$ and $\gamma=1$ , $\mu=0.001$ for the SIRS model. (Left) Phase plane of both solutions, the SIR solution (dashed red) and the SIRS solution (continuous blue), along with the limiting line $S+I=1$ (dashed black). (Right) Comparison of (Top) the solution for the susceptible compartment $S$ of the SIR model (dashed red) against the one of the SIRS model (continuous blue) and (Bottom) the solution of the infectious compartment $I$ of the SIR model (dashed red) against the one of the SIRS model (continuous blue).

4.2.2 Distinguishing between SIR and SIRS models

Given the previous study, it is intuitive to think that, in a practical viewpoint, it may be difficult distinguishing an SIR model from an SIRS model with very small $\mu$ if we look at short times. However, it is theoretically possible. We present here two approaches to do this.

•

Approach 1: One can check analogously to the SIRS case that $\{y,y^{2},\dot{y},y\dot{y}\}$ are linearly independent for the SIR case in any $[a,b)\subset\mathcal{I}$ in its positively invariant set $\tilde{\Omega}=\{(\xi_{1},\xi_{2})^{\mathrm{T}}\in(0,1)^{2}\ :\ \xi_{1}+\xi_{2}\leq 1\}$ . Then, given some observations $y_{(x_{0},\theta_{0})}$ , we can perform Step 3 of Algorithm 3: we know there exist different time instants $t_{1},t_{2},t_{3},t_{4}\in[0,\bar{t}]$ , where $\bar{t}$ is the one in (25), such that there exists a unique solution $\sigma$ fulfilling

\left(\begin{array}[]{cccc}y(t_{1})&y^{2}(t_{1})&\dot{y}(t_{1})&y\dot{y}(t_{1})\\ \\ y(t_{2})&y^{2}(t_{2})&\dot{y}(t_{2})&y\dot{y}(t_{2})\\ \\ y(t_{3})&y^{2}(t_{3})&\dot{y}(t_{3})&y\dot{y}(t_{3})\\ \\ y(t_{4})&y^{2}(t_{4})&\dot{y}(t_{4})&y\dot{y}(t_{4})\end{array}\right)\left(\begin{array}[]{c}\sigma_{1}\\ \sigma_{2}\\ \sigma_{3}\\ \sigma_{4}\end{array}\right)=\left(\begin{array}[]{c}\dfrac{\dot{y}^{2}(t_{1})}{y(t_{1})}-\ddot{y}(t_{1})\\ \dfrac{\dot{y}^{2}(t_{2})}{y(t_{2})}-\ddot{y}(t_{2})\\ \dfrac{\dot{y}^{2}(t_{3})}{y(t_{3})}-\ddot{y}(t_{3})\\ \dfrac{\dot{y}^{2}(t_{4})}{y(t_{4})}-\ddot{y}(t_{4})\end{array}\right).

Given $\sigma$ , we obtain the following conclusions:

–

if $\sigma_{1}=\sigma_{3}=0$ , we confirm that our observations correspond to an SIR model;
–

if $\sigma_{3}\neq 0$ , then we are at an SIRS model;
–

if $\sigma_{1}\neq 0$ and $\sigma_{3}=0$ , we conclude our observations do not match any of the two models.

•

Approach 2: Another way to determine if there is loss of immunity or not consists in considering the previous equality (23), i.e.,

\ddot{y}-\dfrac{\dot{y}}{y}=-\dfrac{\beta\gamma}{k}y^{2}-\dfrac{\beta}{k}y\dot{y},\quad\mbox{which may be rewritten as}\quad\dfrac{\mathrm{d}}{\mathrm{d}t}{\dfrac{\dot{y}}{y}}+\dfrac{\beta\gamma}{k}y+\dfrac{\beta}{k}\dot{y}=0,

i.e., $\left\{\dfrac{\mathrm{d}}{\mathrm{d}t}{\dfrac{\dot{y}}{y}},y,\dot{y}\right\}$ are linearly dependent for the SIR model whenever $I_{\mathrm{SIR}}\not\equiv 0$ . However, this is not true for the SIRS model. Indeed, consider System (12). Let $a_{1},a_{2},a_{3}\in\mathbb{R}$ such that

a_{1}\dfrac{\mathrm{d}}{\mathrm{d}t}{\dfrac{\dot{y}}{y}}+a_{2}y+a_{3}\dot{y}=0.

This is equivalent to

a_{1}\beta(-\beta SI+\mu(1-S-I))+a_{2}kI+a_{3}kI(\beta S-\gamma)=0,

which again reduces to study the linear independence of $1,S,I,SI$ , which we already know that are linearly independent in any $[a,b)\subset\mathcal{I}$ in each set $\Omega_{\theta}$ , $\theta\in\Theta$ , defined in Section 4.1. Therefore, given $y$ , $\dot{y}$ and $\ddot{y}$ , studying the linear dependence of $y,\dot{y},\dfrac{\mathrm{d}}{\mathrm{d}t}\dfrac{\dot{y}}{y}$ may help us determine the model.

5 Some other applications

In this section, we present a series of additional examples to demonstrate the broader applicability of our approach to other epidemiological models. In each example, we illustrate how joined observability-identifiability can be established using our proposed framework, under different modeling assumptions and data limitations. For clarity, we systematically omit the equation of the Recovered compartment in cases where the total population is constant, as it does not affect the theoretical results. To keep the focus on the theoretical aspects, detailed computational steps are omitted but can be derived straightforwardly following the procedures outlined in previous sections.

The SIR model with demography

In [18], the authors show that the following model is neither identifiable nor observable:

\left\{\begin{array}[]{ccl}\dot{S}&=&\delta N-\beta\dfrac{SI}{N}-\delta S,\\[6.00006pt] \dot{I}&=&\beta\dfrac{SI}{N}-(\gamma+\delta)I,\\[6.00006pt] y&=&kI,\end{array}\right.\quad\forall\,t\geq 0,

where $\delta>0$ denotes the birth and death rates, assumed equal, $N$ is the total population (constant), and $\beta$ , $\gamma$ , and $k$ represent the same parameters as in the previous models.

However, this model becomes jointly observable-identifiable when the total population $N$ is known, allowing us to normalize the population so that $S$ and $I$ now represent fractions of the total population:

\left\{\begin{array}[]{ccl}\dot{S}&=&\delta-\beta SI-\delta S,\\ \dot{I}&=&\beta SI-(\gamma+\delta)I,\\ y&=&kI,\end{array}\right.\quad\forall\,t\geq 0,

(26)

for some initial condition $(S(0),I(0))^{\mathrm{T}}\in[0,1]^{2}$ .

For this normalized system, we set $\Omega=\{(\xi_{1},\xi_{2})^{\mathrm{T}}\in[0,1]^{2}\,:\,\xi_{1}+\xi_{2}\leq 1\}$ and $\Theta:=(0,1]\times(0,+\infty)^{3}$ for $\theta=(k,\beta,\gamma,\delta)$ . Performing computations analogous to those in Section 4.1 during Step 1 and Step 2, we obtain the following relations:

y=kI,\quad\dot{y}=(\beta S-(\gamma+\delta))y\quad\implies\quad S=\dfrac{1}{\beta}\left(\dfrac{\dot{y}}{y}+\gamma+\delta\right),\quad I=\dfrac{y}{k},

assuming $y\not\equiv 0$ , and

\ddot{y}-\dfrac{\dot{y}^{2}}{y}=y\delta(\beta-\delta-\gamma)-y^{2}\dfrac{\beta}{k}(\gamma+\delta)-\dot{y}\delta-y\dot{y}\dfrac{\beta}{k}.

Then, we consider $d_{1}=1$ , $d_{1}^{\prime}=2$ , $p=1$ , $q_{1}=q=4$ ,

g_{0}(y,\dot{y},\ddot{y})=\ddot{y}-\dfrac{\dot{y}^{2}}{y},\ g_{1}(y,\dot{y},\ddot{y})=y,\ g_{2}(y,\dot{y},\ddot{y})=-y^{2},\ g_{3}(y,\dot{y},\ddot{y})=-\dot{y},\ g_{4}(y,\dot{y},\ddot{y})=-y\dot{y}

and

r(\theta)=\left(\delta(\beta-\gamma-\delta),\dfrac{\beta}{k}(\gamma+\delta),\delta,\dfrac{\beta}{k}\right).

With these expressions, we confirm that the assumptions in Step 1 and Step 2 of Algorithm 3 are satisfied, taking parameters in $\Theta$ and initial conditions consistent with the family $\{\Omega_{\theta}\}_{\theta\in\Theta}$ , with $\Omega_{\theta}:=\{(\xi_{1},\xi_{2})\in[0,1)\times(0,1]\ :\ \xi_{1}+\xi_{2}\leq 1\}\setminus\{P_{\mathrm{e}}(\theta)\}$ , for each $\theta\in\Theta$ , where $P_{\mathrm{e}}(\theta)$ denotes the equilibrium point associated with $\theta$ , and $\Omega_{\theta}$ is positively invariant with respect to the ODE system of (26) when we consider $\theta$ as parameter vector.

Therefore, we conclude that System (26) is jointly observable-identifiable on $\Gamma_{\Theta}$ in any $[a,b)\subset\mathcal{I}$ .

The SIRV model

We now analyze a Susceptible-Infectious-Recovered-Vaccinated (SIRV) model, where individuals gain permanent immunity after recovery, as in the SIR model. Additionally, we assume the existence of a perfect vaccine that provides permanent immunity but is only effective for susceptible individuals. Importantly, this vaccine has no effect on infectious or recovered individuals. However, since we cannot distinguish susceptible individuals from infectious or recovered ones, all compartments are vaccinated indiscriminately. The observable data are the rate of vaccinated individuals (including both effective and ineffective vaccinations). Hence, we consider the following system, where $S$ , $I$ , and $V$ are fractions of the total population, assumed to be known:

\left\{\begin{array}[]{ccl}\dot{S}&=&-\beta SI-\nu S,\\ \dot{I}&=&\beta SI-\gamma I,\\ \dot{V}&=&\nu S,\\ y&=&\nu(1-V),\end{array}\right.\quad\forall\,t\geq 0,

(27)

for some initial condition $(S(0),I(0),V(0))^{\mathrm{T}}\in[0,1]^{3}$ . We assume that the initial condition and all parameters $\beta$ , $\gamma$ , and $\nu$ are unknown, where $\beta$ and $\gamma$ are as defined in previous models, and $\nu>0$ (days^-1) represents the vaccination rate.

For this system, we define $\Omega=\{(\xi_{1},\xi_{2},\xi_{3})^{\mathrm{T}}\in[0,1]^{3}\ :\ \xi_{1}+\xi_{2}+\xi_{3}\leq 1\}$ and $\Theta:=(0,+\infty)^{3}$ for $\theta=(\beta,\gamma,\nu)$ . After some computations, we obtain the following relations:

y=\nu(1-V),\quad\dot{y}=-\nu^{2}S,\quad\ddot{y}=-\dot{y}(\beta I+\nu)\quad\implies\quad S=-\dfrac{\dot{y}}{\nu^{2}},\quad I=-\dfrac{1}{\beta}\left(\dfrac{\ddot{y}}{\dot{y}}+\nu\right),\quad V=1-\dfrac{y}{\nu},

assuming $\dot{y}\not\equiv 0$ , and

y^{(3)}-\dfrac{\ddot{y}^{2}}{\dot{y}}=-\dot{y}\nu\gamma-\dot{y}^{2}\dfrac{\beta}{\nu}-\ddot{y}\gamma-\dot{y}\ddot{y}\dfrac{\beta}{\nu^{2}}.

Then, we consider $d_{1}=2$ , $d_{1}^{\prime}=3$ , $p=1$ , $q_{1}=q=4$ ,

g_{0}(y,\dot{y},\ddot{y},y^{(3)})=y^{(3)}-\dfrac{\ddot{y}^{2}}{\dot{y}},\quad g_{1}(y,\dot{y},\ddot{y},y^{(3)})=-\dot{y},\quad g_{2}(y,\dot{y},\ddot{y},y^{(3)})=-\dot{y}^{2},

g_{3}(y,\dot{y},\ddot{y},y^{(3)})=-\ddot{y},\quad g_{4}(y,\dot{y},\ddot{y},y^{(3)})=-\dot{y}\ddot{y},

and

r(\theta)=\left(\nu\gamma,\dfrac{\beta}{\nu},\gamma,\dfrac{\beta}{\nu^{2}}\right).

Following the procedures established in earlier sections, we confirm that the assumptions of Step 1 and Step 2 in Algorithm 3 are met, with parameters in $\Theta$ and initial conditions in the revised set $\Omega=\{(\xi_{1},\xi_{2},\xi_{3})\in(0,1]\times(0,1)\times[0,1)\ :\ \xi_{1}+\xi_{2}+\xi_{3}\leq 1\}$ , where $\Omega$ does not depend on $\theta\in\Theta$ and is positively invariant under the ODE system defined by (27).

Thus, we conclude that System (27) is jointly observable-identifiable on $\Omega\times\Theta$ in any $[a,b)\subset\mathcal{I}$ .

The SIR model (with a different output)

As discussed in Section 4.2, the SIR model observed through a fraction of infected individuals (i.e., System (21)) is observable and identifiable, but not jointly observable-identifiable. Here, we consider the same SIR model but with a different observation: the instantaneous incidence rate. Specifically, we examine the following system:

\left\{\begin{array}[]{lcl}\dot{S}&=&-\beta SI,\\ \dot{I}&=&\beta SI-\gamma I,\\ y&=&\beta SI,\end{array}\right.

(28)

for some initial condition $(S(0),I(0))\in\Omega=\{(\xi_{1},\xi_{2})^{\mathrm{T}}\in[0,1]^{2}\ :\ \xi_{1}+\xi_{2}\leq 1\}$ , which is positively invariant with respect to the system of ODEs given in (28). In this scenario, we assume that both $\beta$ and $\gamma$ are unknown parameters, defining $\Theta:=(0,\infty)^{2}$ for $\theta=(\beta,\gamma)$ . Following an approach similar to previous cases, we calculate the first derivative of $y$ to express $S$ and $I$ in terms of $y$ and $\dot{y}$ :

\left\{\begin{array}[]{lcl}y&=&\beta SI,\\ \dot{y}&=&-y\beta I+y\beta S-\gamma y.\end{array}\right.

Notice that this system is nonlinear in $S$ and $I$ , making it challenging to directly apply Step 1 of Algorithm 3. To address this, let us differentiate once more:

\ddot{y}=-\dot{y}\beta I-2y^{2}\beta+y\beta\gamma I+\dot{y}\beta S-\gamma\dot{y},

which results in a system that is linear in $S$ and $I$ . Therefore, we obtain:

\left\{\begin{array}[]{lcl}\dot{y}&=&(-\beta I+\beta S-\gamma)y,\\ \ddot{y}&=&(-\beta I+\beta S-\gamma)\dot{y}-2y^{2}\beta+y\beta\gamma I,\end{array}\right.\implies\left\{\begin{array}[]{lcl}S&=&\dfrac{1}{\beta}\dfrac{\dot{y}}{y}+\dfrac{1}{\beta\gamma}\left(\dfrac{\ddot{y}}{y}-\dfrac{\dot{y}^{2}}{y^{2}}\right)+\dfrac{2}{\gamma}y+\dfrac{\gamma}{\beta},\\[10.00002pt] I&=&\dfrac{1}{\beta\gamma}\left(\dfrac{\ddot{y}}{y}-\dfrac{\dot{y}^{2}}{y^{2}}\right)+\dfrac{2}{\gamma}y.\end{array}\right.

Assuming $y\not\equiv 0$ , we avoid further differentiation by substituting the expressions of $S$ and $I$ into $y=\beta SI$ , yielding the following equation after simplification:

-\left(\dfrac{\mathrm{d}}{\mathrm{d}t}\dfrac{\dot{y}}{y}\right)^{2}=\gamma\dfrac{\dot{y}}{y}\dfrac{\mathrm{d}}{\mathrm{d}t}\dfrac{\dot{y}}{y}+2\beta\gamma\dot{y}+\beta\gamma^{2}y+4\beta y\dfrac{\mathrm{d}}{\mathrm{d}t}\dfrac{\dot{y}}{y}+4\beta^{2}y^{2}+\gamma^{2}\dfrac{\mathrm{d}}{\mathrm{d}t}\dfrac{\dot{y}}{y}.

Then, we set $d_{1}=d_{1}^{\prime}=2$ , $p=1$ , $q_{1}=q=6$ ,

g_{0}(y,\dot{y},\ddot{y})=-\left(\dfrac{\mathrm{d}}{\mathrm{d}t}\dfrac{\dot{y}}{y}\right)^{2},\quad g_{1}(y,\dot{y},\ddot{y})=\dfrac{\dot{y}}{y}\dfrac{\mathrm{d}}{\mathrm{d}t}\dfrac{\dot{y}}{y},\quad g_{2}(y,\dot{y},\ddot{y})=\dot{y},\quad g_{3}(y,\dot{y},\ddot{y})=y,

g_{4}(y,\dot{y},\ddot{y})=y\dfrac{\mathrm{d}}{\mathrm{d}t}\dfrac{\dot{y}}{y},\quad g_{5}(y,\dot{y},\ddot{y})=y^{2},\quad g_{6}(y,\dot{y},\ddot{y})=\dfrac{\mathrm{d}}{\mathrm{d}t}\dfrac{\dot{y}}{y},

and

r(\theta)=\left(\gamma,2\beta\gamma,\beta\gamma^{2},4\beta,4\beta^{2},\gamma^{2}\right).

We confirm that the assumptions of Step 1 and Step 2 in Algorithm 3 hold, with parameters in $\Theta$ and initial conditions in $\Omega=\{(\xi_{1},\xi_{2})\in(0,1)^{2}\ :\ \xi_{1}+\xi_{2}\leq 1\}$ , where $\Omega$ does not depend on $\theta\in\Theta$ and is positively invariant with respect to the system of ODEs of System (28).

Consequently, we conclude that System (28) is jointly observable-identifiable on $\Omega\times\Theta$ in any $[a,b)\subset\mathcal{I}$ .

The SIV model with demography and two outputs

Here, we consider a Susceptible-Infectious-Vaccinated (SIV) model with demographic dynamics. We assume a perfect vaccine that is effective only for susceptible individuals, though infectious individuals are also vaccinated indiscriminately, as in System (27). In this case, however, we model a permanent infection, which can serve as a simplified representation of long-term infections, such as those associated with certain sexually transmitted diseases (e.g., human papillomavirus (HPV)). We assume two observable quantities: the rate of all vaccinated individuals (both effective and ineffective) and the rate of natural deaths. These deaths are attributed to natural population dynamics and are assumed to be routinely monitored by authorities. Hence, we consider the following system, where $S$ , $I$ , and $V$ represent the fraction of susceptible, infectious, and vaccinated individuals, respectively:

\left\{\begin{array}[]{lcl}\dot{S}&=&A-\beta SI-(\nu+\delta)S,\\ \dot{I}&=&\beta SI-\delta I,\\ \dot{V}&=&\nu S-\delta V,\\ y&=&\left(\begin{array}[]{c}y_{1}\\ y_{2}\end{array}\right)=\left(\begin{array}[]{c}\nu(S+I)\\ \delta(S+I+V)\end{array}\right),\end{array}\right.\quad\forall\,t\geq 0,

(29)

for some initial condition $(S(0),I(0),V(0))\in[0,\infty)^{3}$ . All the initial condition and the parameters $A$ , $\beta$ , $\delta$ , and $\nu$ are assumed to be unknown. Here, $A>0$ represents a constant recruitment rate, $\delta>0$ is the death rate, and $\beta$ and $\nu$ have the same definitions as in System (27).

The natural positively invariant set for this system is $\Omega=\{(\xi_{1},\xi_{2},\xi_{3})^{\mathrm{T}}\in\left[0,A/\delta\right]^{3}\ :\ \xi_{1}+\xi_{2}+\xi_{3}\leq A/\delta\}$ , and $\Theta:=(0,\infty)^{4}$ , for $\theta=(A,\beta,\delta,\nu)$ .

We first observe that this system is not observable (and hence not jointly observable-identifiable) if only $y_{1}$ is used as the observation, similar to System (27). Since $y_{1}$ depends solely on $S$ and $I$ , we cannot uniquely determine $V_{0}$ . Therefore, an additional observation is required. By including the rate of natural deaths, $y_{2}=\delta(S+I+V)$ , we leverage data that should be routinely available, making it reasonable to assume accessibility to this information.

Similarly, relying on $y_{2}$ alone is insufficient, as $\dot{y}_{2}=A\delta-\delta y_{2}$ yields no further parameter information through differentiation, unlike previous models.

After some computations, differentiating $y_{1}$ twice and $y_{2}$ once, we obtain the following expressions:

\left\{\begin{array}[]{lcl}y_{1}&=&\nu(S+I),\\ \dot{y}_{1}&=&\nu A-\nu^{2}S-\delta y_{1},\\ y_{2}&=&\delta(S+I+V)\end{array}\right.\implies\left\{\begin{array}[]{lcl}S&=&\dfrac{A}{\nu}-\dfrac{\delta}{\nu^{2}}y_{1}-\dfrac{\dot{y}_{1}}{\nu^{2}},\\[10.00002pt] I&=&\left(\dfrac{\nu+\delta}{\nu^{2}}\right)y_{1}-\dfrac{A}{\nu}+\dfrac{\dot{y}_{1}}{\nu^{2}},\\[10.00002pt] V&=&\dfrac{y_{2}}{\delta}-\dfrac{y_{1}}{\nu},\end{array}\right.

and

\left\{\begin{array}[]{lcl}\dot{y}_{1}^{2}&=&A\nu^{2}\dfrac{\delta\nu-A\beta}{\beta}+\nu\left(A\nu+2A\delta-\dfrac{\delta\nu^{2}+\delta^{2}\nu}{\beta}\right)y_{1}+\nu\left(2A-\dfrac{\nu^{2}+2\delta\nu}{\beta}\right)\dot{y}_{1}\\ &&-\delta(\nu+\delta)y_{1}^{2}-(\nu+2\delta)y_{1}\dot{y}_{1}-\dfrac{\nu^{2}}{\beta}\ddot{y}_{1},\\ \dot{y}_{2}&=&A\delta-\delta y_{2}.\end{array}\right.

Focusing on the equation $\dot{y}_{2}=A\delta-\delta y_{2}$ , we observe that $\{1,y_{2}\}$ , when $y_{2}$ is non-constant, allows us to determine $A$ and $\delta$ by solving the system:

\begin{pmatrix}1&y_{2}(t_{1})\\ 1&y_{2}(t_{2})\end{pmatrix}\begin{pmatrix}\sigma_{1}\\ \sigma_{2}\end{pmatrix}=\begin{pmatrix}\dot{y}_{2}(t_{1})\\ \dot{y}_{2}(t_{2})\end{pmatrix},

such that $\delta=-\sigma_{2}$ and $A=-\sigma_{1}/\sigma_{2}$ . Substituting, we can rewrite:

\left\{\begin{array}[]{rcl}\dot{y}_{1}^{2}+\sigma_{2}^{2}y_{1}^{2}-2y_{1}\dot{y}_{1}\sigma_{2}&=&A\nu^{2}\dfrac{\delta\nu-A\beta}{\beta}+\nu\left(A\nu+2A\delta-\dfrac{\delta\nu^{2}+\delta^{2}\nu}{\beta}\right)y_{1}\\[10.00002pt] &&+\nu\left(2A-\dfrac{\nu^{2}+2\delta\nu}{\beta}\right)\dot{y}_{1}-\nu y_{1}\left(-\sigma_{2}y_{1}+\dot{y}_{1}\right)-\dfrac{\nu^{2}}{\beta}\ddot{y}_{1},\\ \dot{y}_{2}&=&A\delta-\delta y_{2}.\end{array}\right.

This reformulation ensures the linear independence condition in Step 2 of Algorithm 3. Then, we consider $d_{1}=1$ , $d_{2}=0$ , $d_{1}^{\prime}=2$ , $d_{2}^{\prime}=1$ , $p=2$ , $q_{1}=5$ , $q_{2}=2$ ,

$\displaystyle g_{1,0}(y_{1},\dot{y}_{1},\ddot{y}_{1},y_{2},\dot{y}_{2})$	$\displaystyle=$	$\displaystyle\dot{y}_{1}^{2}+\sigma_{2}y_{1}^{2}-2y_{1}\dot{y}_{1}\sigma_{2},\quad g_{1,1}(y_{1},\dot{y}_{1},\ddot{y}_{1},y_{2},\dot{y}_{2})=1,$
$\displaystyle g_{1,2}(y_{1},\dot{y}_{1},\ddot{y}_{1},y_{2},\dot{y}_{2})$	$\displaystyle=$	$\displaystyle y_{1},\quad g_{1,3}(y_{1},\dot{y}_{1},\ddot{y}_{1},y_{2},\dot{y}_{2})=\dot{y}_{1},$
$\displaystyle g_{1,4}(y_{1},\dot{y}_{1},\ddot{y}_{1},y_{2},\dot{y}_{2})$	$\displaystyle=$	$\displaystyle-y_{1}\left(-\sigma_{2}y_{1}+\dot{y}_{1}\right),\quad g_{1,5}(y_{1},\dot{y}_{1},\ddot{y}_{1},y_{2},\dot{y}_{2})=-\ddot{y}_{1},$
$\displaystyle g_{2,0}(y_{1},\dot{y}_{1},\ddot{y}_{1},y_{2},\dot{y}_{2})$	$\displaystyle=$	$\displaystyle\dot{y}_{2},\quad g_{2,1}(y_{1},\dot{y}_{1},\ddot{y}_{1},y_{2},\dot{y}_{2})=1,\quad g_{2,2}(y_{1},\dot{y}_{1},\ddot{y}_{1},y_{2},\dot{y}_{2})=-y_{2},$

and

r(\theta_{0})=\left(A\nu^{2}\dfrac{\delta\nu-A\beta}{\beta},\nu\left(A\nu+2A\delta-\dfrac{\delta\nu^{2}+\delta^{2}\nu}{\beta}\right),\nu\left(2A-\dfrac{\nu^{2}+2\delta\nu}{\beta}\right),\nu,\dfrac{\nu^{2}}{\beta},A\delta,\delta\right).

Given this, the assumptions in Step 1 and Step 2 of Algorithm 3 hold for parameters in $\Theta$ and initial conditions consistent with the family $\{\Omega_{\theta}\}_{\theta\in\Theta}$ , such that, for $\theta\in\Theta$ , $\Omega_{\theta}:=\{(\xi_{1},\xi_{2},\xi_{3})\in[0,A/\delta)\times(0,A/\delta)\times[0,A/\delta)\ :\ \xi_{1}+\xi_{2}+\xi_{3}<A/\delta\}\setminus\{P_{\mathrm{e}}(\theta)\}$ , where $P_{\mathrm{e}}(\theta)$ is the equilibrium point of the system and $\Omega_{\theta}$ is positively invariant under the ODEs of System (29) when considering $\theta$ as the parameter vector.

Therefore, we conclude that System (29) is jointly observable-identifiable on $\Gamma_{\Theta}$ in any $[a,b)\subset\mathcal{I}$ .

6 Numerical illustration

In this section, we present numerical experiments designed to illustrate the application of the theoretical framework developed in Sections 3 and 4 for observability, identifiability and joined observability-identifiability in epidemiological models. This numerical analysis also aims to implement the theoretical methods in some cases and demonstrate their practical feasibility under varying conditions. Specifically, we analyze two examples: one based on an SIRS model and another on an SIR model, each chosen to exhibit distinct epidemiological dynamics:

•

Case 1: We consider an SIRS model in which the endemic equilibrium is globally asymptotically stable (except for the invariant manifold $I\equiv 0$ , as it is stable for the disease-free equilibrium). The solution oscillates when converging to the endemic equilibrium, as in the example presented in Figure 1. We set the following parameters and initial conditions: $k_{0}=0.3,\ \beta_{0}=0.25,\ \gamma_{0}=0.1,\ \mu_{0}=0.05,\ S_{0}=0.9,$ and $I_{0}=0.1$ . Thus, the basic reproduction number is $\mathcal{R}_{0}=2.5$ .
•

Case 2: We consider an SIR model that exhibits an epidemic peak. The chosen parameters and initial conditions are $k_{0}=0.3,\ \beta_{0}=0.25,\ \gamma_{0}=0.1,\ S_{0}=0.9,$ and $I_{0}=0.1$ . Here as well, $\mathcal{R}_{0}=2.5$ .

To simulate these cases, we approximate each solution using a fourth-order Runge-Kutta algorithm and synthetically generate observations, $y=k_{0}I$ , at each time-step. In this initial study, we assume noise-free data to establish a clear baseline for comparison. We then analyze two scenarios:

•

Scenario 1: We assume that continuous observations of $y$ are available over the interval $[0,T_{\max}]$ , with exact knowledge or computability of its derivatives (enabled by the analyticity of the system). We apply the recovery procedure outlined in Algorithm 4 to determine the original parameters and initial condition, utilizing the full interval $[0,T_{\max}]$ .
•

Scenario 2: In realistic situations, data are typically discrete, often recorded at daily intervals by public health authorities (e.g., [68] and [69]). To mimic this, we extract simulated data once per day (at the same time each day), referred to here as daily data, and apply the Ordinary Least Squares method to estimate the unknown parameters.

Scenario 1 provides a controlled environment for testing recovery procedures and Scenario 2 reflects more realistic, discrete data conditions. For practical purposes, we will examine Case 1 under Scenario 1 and Case 2 under Scenario 2. Preliminary results indicated similar outcomes when interchanging cases between scenarios, allowing us to streamline the analysis. In Scenario 1, we focus on the procedure in Algorithm 4, omitting that of Algorithm 3 due to the former’s superior performance.

To integrate numerically the ODE systems, we utilize a fourth-order Runge-Kutta method with a time-step of $h=2^{-5}$ days, extending up to a maximum simulation time of $T_{\max}=5$ days. Although the maximum time may seem low, it was sufficient to capture the dynamics required for convergence in Scenario 2. Furthermore, we employ the extended SIRS model from Section 4.2, assuming no prior knowledge regarding whether $\mu_{0}\neq 0$ or $\mu_{0}=0$ .

6.1 Scenario 1: Linear systems for continuous observations

In this section, we perform numerical tests for Case 1, assuming continuous observation of $y=k_{0}I$ over $[0,T_{\max}]$ , i.e., we observe $y$ at every time point within this interval, and assume that we know or can compute its exact successive derivatives.

To implement the numerical tests, we first need to construct the necessary data. From an implementation perspective, we obtain the values of $y$ and its derivatives at each time-step of the numerical scheme, specifically at each point in the time vector $t_{\mathrm{vec}}=(t_{i})_{i\in\{0,\dots,N\}}$ , where $t_{i}=ih$ and $N=T_{\max}/h\in\mathbb{N}$ . The data generation steps are as follows:

1.

We approximate the solutions $S$ and $I$ of the SIR or SIRS system with a small time-step $h$ using the fourth-order Runge-Kutta scheme, and define our observations as $y=k_{0}I$ .
2.

We approximate the first derivative of $y$ by $\dot{y}=k_{0}\dot{I}=k_{0}I(\beta_{0}S-\gamma_{0})$ .
3.

For higher-order derivatives, we use the following linear equation in terms of parameters $\sigma_{i}$ , $i\in\{1,2,3,4\}$ :

$\ddot{y}=\dfrac{\dot{y}^{2}}{y}-\sigma_{1}y-\sigma_{2}y^{2}-\sigma_{3}\dot{y}-\sigma_{4}y\dot{y},$

where $\sigma_{1}=\mu_{0}(\gamma_{0}-\beta_{0})$ , $\sigma_{2}=\beta_{0}(\gamma_{0}+\mu_{0})/k_{0}$ , $\sigma_{3}=\mu_{0}$ , and $\sigma_{4}=\beta_{0}/k_{0}$ . Using this relationship, we iteratively approximate higher-order derivatives of $y$ through lower-order derivatives.

With these generated data, we apply the method in Algorithm 4 to estimate the system parameters and initial condition. This method requires selecting a time $\tilde{t}\in t_{\mathrm{vec}}$ such that the following system has a unique solution (see (10)-(11)):

\left(\begin{array}[]{cccc}\dfrac{\mathrm{d}^{k}y}{\mathrm{d}t^{k}}&\dfrac{\mathrm{d}^{k}y^{2}}{\mathrm{d}t^{k}}&\dfrac{\mathrm{d}^{k}\dot{y}}{\mathrm{d}t^{k}}&\dfrac{\mathrm{d}^{k}y\dot{y}}{\mathrm{d}t^{k}}\end{array}\right)_{\begin{subarray}{l}t=\tilde{t}\\ k=0,1,2,3\end{subarray}}\sigma^{\mathrm{T}}=\left(\begin{array}[]{c}\dfrac{\mathrm{d}^{k}}{\mathrm{d}t^{k}}\left(\dfrac{\dot{y}^{2}}{y}-\ddot{y}\right)\end{array}\right)_{\begin{subarray}{l}t=\tilde{t}\\ k=0,1,2,3\end{subarray}}.

(30)

Instead of directly estimating $(k_{0},\beta_{0},\gamma_{0},\mu_{0})$ , we focus on estimating

\sigma=(\sigma_{1},\sigma_{2},\sigma_{3},\sigma_{4})=\left(\mu_{0}(\gamma_{0}-\beta_{0}),\dfrac{\beta_{0}}{k_{0}}(\gamma_{0}+\mu_{0}),\mu_{0},\dfrac{\beta_{0}}{k_{0}}\right).

The initial conditions $S_{0}$ and $I_{0}$ are then computed as

S_{0}=\dfrac{1}{\beta_{0}}\left(\dfrac{\dot{y}(0)}{y(0)}+\gamma_{0}\right)\quad\text{and}\quad I_{0}=\dfrac{y(0)}{k_{0}}.

Knowing the observation at time 0 avoids the computational challenges of backward integration.

Assuming known values of $y(0)=0.03$ and $\dot{y}(0)=0.00375$ , we aim to calibrate

(\sigma_{1},\sigma_{2},\sigma_{3},\sigma_{4})=(-0.0075,\ 0.125,\ 0.05,\ 0.8\overset{\_}{3})\quad\text{and}\quad(S_{0},I_{0})=(0.9,0.1).

To evaluate the method’s performance, we conducted 161 experiments, corresponding to $N+1=T_{\max}/h+1=161$ . For each time point in $t_{\mathrm{vec}}$ , we examined:

•

The relative error between each computed component of $\sigma$ and its exact value. These errors were consistently small, with a maximum of the order of $10^{-13}$ , confirming the high accuracy of the method.
•

The determinant of the matrix in (30), which should be non-zero. The values obtained were of the order of $10^{-17}$ .
•

The condition number of the matrix in (30) (based on the $L^{2}$ -norm). All the obtained values range between $10^{3}$ and $5\times 10^{5}$ .
•

The computational time (in seconds). Performing all tests took approximately $10^{-2}$ seconds.

These results indicate that the method provides accurate parameter estimates despite the small determinant values and the matrix conditioning. The computational efficiency demonstrated, completing 161 tests in minimal time, makes this approach highly suitable for real-time applications.

Similar results were observed for Case 2, which are omitted here for brevity.

6.2 Scenario 2: OLS method for daily observations

In this section, we analyze Case 2 under the assumption of daily, noise-free observations from the deterministic model. Daily observations reflect realistic data collection practices in epidemiology, such as those by public health authorities. Given the discrete nature of the data, calculating derivatives directly can result in inaccuracy (if differentiated numerically) or bias (if interpolated prior to differentiation). In this context, we employ the Ordinary Least Squares (OLS) method to estimate the values of $k_{0}$ , $\beta_{0}$ , $\gamma_{0}$ , $\mu_{0}$ , and $S_{0}$ . Let $\theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4},\theta_{5})$ represent the parameter vector, where we calibrate 4 parameters ( $k_{0}$ , $\beta_{0}$ , $\gamma_{0}$ , and $\mu_{0}$ ) and 1 initial condition ( $S_{0}$ ). The goal is to find $\theta_{\mathrm{OLS}}=(k_{\mathrm{OLS}},\beta_{\mathrm{OLS}},\gamma_{\mathrm{OLS}},\mu_{\mathrm{OLS}},S_{0,\mathrm{OLS}})$ by solving

\theta_{\mathrm{OLS}}=\underset{\theta\in{\Theta}}{\mathrm{argmin}}\sum_{i=0}^{T_{\max}}A\left(\theta_{1}I(t_{i};\theta)-y(t_{i})\right)^{2},\quad t_{i}=i,

where $I(t;\theta)$ is the solution to the following system of ODEs:

\left\{\begin{array}[]{lcl}\dot{S}(t;\theta)&=&-\theta_{2}S(t;\theta)I(t;\theta)+\theta_{4}(1-S(t;\theta)-I(t;\theta)),\\[6.00006pt] \dot{I}(t;\theta)&=&\theta_{2}S(t;\theta)I(t;\theta)-\theta_{3}I(t;\theta),\end{array}\right.

with $(S(0;\theta),I(0;\theta))=\left(\theta_{5},\dfrac{y(0)}{\theta_{1}}\right)$ . Here, we define a new feasible set $\Theta^{*}=\Theta\times[0,1)$ . To enhance sensitivity to small deviations, a scaling factor $A=10^{14}$ is introduced in the objective function.

The choice of $T_{\max}=5$ days (using daily data) was based on preliminary testing, where this time-frame was found sufficient for reliable parameter estimation despite its brevity.

For this experiment, we use the MATLAB function lsqcurvefit to solve the OLS problem, adjusting settings for increased accuracy (see [50]). Specifically, we set the maximum evaluations to $10^{6}$ , iterations to $5\times 10^{5}$ , step tolerance to $10^{-15}$ , and function tolerance to $10^{-17}$ .

We impose the following bounds for the parameters, with subscripts $m$ and $M$ indicating the minimum and maximum values, respectively: $[k_{m},k_{M}]=[\min_{t}\{y\},1]$ , $[\beta_{m},\beta_{M}]=[10^{-2},3]$ , $[\gamma_{m},\gamma_{M}]=[10^{-2},1]$ , $[\mu_{m},\mu_{M}]=[0,1]$ , and $[S_{m},S_{M}]=[0,1-10^{-10}]$ . The bounds for $\beta$ are informed by studies on $\mathcal{R}_{0}$ estimates for diseases such as smallpox, pertussis, or COVID-19 (e.g., [23], [43], [45], [46], [66], [67]).

We initialize the OLS algorithm several times with random initial conditions, denoted as I.C., drawn from a uniform distribution $\mathcal{U}(0,1)$ . Each initial condition is generated as:

\text{I.C.}=(k_{m},\beta_{m},\gamma_{m},\mu_{m},S_{m})^{\mathrm{T}}+(\rho_{1}(k_{M}-k_{m}),\rho_{2}(\beta_{M}-\beta_{m}),\rho_{3}(\gamma_{M}-\gamma_{m}),\rho_{4}(\mu_{M}-\mu_{m}),\rho_{5}(S_{M}-s_{m}))^{\mathrm{T}},

where $\rho_{i}\sim\mathcal{U}(0,1)$ . For conciseness, we report results from 5 representative tests. Table 2 includes the following outcomes: Test (test number), I.C. (initial condition for lsqcurvefit), Abs. error $\theta_{0}$ (absolute error between the computed solution and the true parameter vector $\theta_{0}$ ), Obj. value (final objective function value), and Time (s) (computational time for each test in seconds).

Test	I.C.	Abs. error $\theta_{0}$	Obj. value	Time (s)
1	(0.545 1.967 0.414 0.82 0.718)	(0.019 0.015 6.8e-6 4.7e-6 0.052)	6.324e-12	116.517
2	(0.97 1.599 0.332 0.106 0.611)	(0.226 0.189 8.9e-12 9.3e-13 0.387)	7.657e-19	1.904
3	(0.785 1.276 0.1 0.266 0.154)	(0.25 0.208 4e-11 3.6e-12 0.409)	9.036e-18	1.546
4	(0.686 0.874 0.675 0.695 0.068)	(0.7 0.286 0.270 0.096 0.024)	0.003	621.808
5	(0.277 0.68 0.671 0.844 0.344)	(0.7 0.287 0.271 0.097 0.024)	0.006	1380.792

Table 1: Results for Scenario 2 applied to Case 2 using the OLS method.

Test	Approx. of $\theta_{0}$	Approx. of $\tilde{\theta}_{0}$
1	(0.319 0.265 0.1 4.7e-6 0.848)	(0.1 0.833 0.225)
2	(0.526 0.439 0.1 9.3e-13 0.513)	(0.1 0.833 0.225)
3	(0.55 0.458 0.1 3.6e-12 0.491)	(0.1 0.833 0.225)
4	(1 0.536 0.370 0.096 0.924)	(0.370 0.536 0.495)
5	(1 0.537 0.371 0.097 0.924)	(0.371 0.537 0.496)

Table 2: Results for Scenario 2 applied to Case 2 using the OLS method: Approximations of the parameters.

The target parameter vector is

\theta_{0}=(k_{0},\beta_{0},\gamma_{0},\mu_{0},S_{0})=(0.3,\ 0.25,\ 0.1,\ 0,\ 0.9),

but as noted in Section 4.2, we can only determine specific combinations: $\gamma_{0}$ , $\beta_{0}/k_{0}$ , $\beta S_{0}$ , $kI_{0}$ . Defining $\tilde{\theta}_{0}=(\gamma_{0},\beta_{0}/k_{0},\beta_{0}S_{0})$ ,

\tilde{\theta}_{0}=(\gamma_{0},\beta_{0}/k_{0},\beta_{0}S_{0})=(0.1,0.8\overset{\_}{3},0.225).

Table 2 compares approximations for $\theta_{0}$ and $\tilde{\theta}_{0}$ , showing accurate convergence for $\gamma_{0}$ and $\mu_{0}$ in Tests 1, 2 and 3. Despite inaccuracies in $k_{0}$ , $\beta_{0}$ , and $S_{0}$ in these tests, the objective function values are low, indicating that the algorithm converged successfully for the key identifiable combinations. Most of the tests not presented here converged similarly.

In Tests 4 and 5, the parameters reach a similar solution, representing an SIRS model with

k\approx 1,\,\beta\approx 0.54,\,\gamma\approx 0.37,\,\mu\approx 0.1,\,\text{and}\,S(0)\approx 0.924.

These results further illustrate the method’s capacity to handle challenging scenarios where early-stage data provide limited discriminatory power. The amplified error factor of $10^{14}$ implies that the non-amplified error in these cases is of the order of $10^{-17}$ , illustrating the difficulty in distinguishing between SIR and SIRS models in early stages when only a portion of infected individuals is observed, as noted in Section 4.2.1. Several tests not presented here reached similar solutions. The rest of the tests showed different results.

For Case 1, similar results were observed, with certain tests producing parameter sets indicative of an SIR model (i.e., with $\mu\approx 0$ ). Nonetheless, these tests consistently preserved the identifiable combinations $\gamma$ , $\beta/k$ and $\beta S_{0}$ , resulting in undistinguishable SIR models. These results are not presented here for brevity.

7 Conclusions

This work addresses the problems of observability, identifiability and joined identifiability-observability in a wide class of systems of ODEs that covers classical epidemiological models, considering idealized continuous, noise-free observations. This problem was previously proposed in [44], where the basis of the approach was developed in a general context. In Section 3, we extended the theoretical framework in [44] considering parameter-dependent sets of initial conditions, and presented Algorithms 3, 3 and 4, offering different constructive methods to compute unknown parameters and/or initial conditions based on this framework. These methods rely on solving a series of linear systems at well chosen time instants, providing a practical approach to recovering model parameters. Additionally, if observations are available at time 0, the initial condition can be determined directly; otherwise, it can be inferred by integrating backwards using the estimated parameters. Furthermore, our approach can be extended to models with piecewise constant parameters by allowing parameter updates at predetermined time instants.

One of the primary ideas underlying our theoretic framework is the establishment of the linear independence of some sets of functions, an approach that has not been much addressed, and the known results are framed in more restricted settings (e.g., rational models, or under the knowledge of data at initial time). This technique, although challenging depending on the sets of functions, has proved to be essential for deriving the joined observability-identifiability results. Moreover, this method reduces the dependency on high-order derivatives, offering simpler computations and more robust algorithms against noise compared to classical techniques such as the Hermann-Krener matrix rank condition.

In Section 4, we applied this framework to an SIRS model (with $\mu>0$ ) when only a fraction of infected individuals is observed, under the assumption of continuous, noise-free observations, and demonstrated its theoretical joined observability-identifiability. Specifically, we showed that observing only a fraction of the infected population is sufficient to determine all model parameters ( $k$ , $\beta$ , $\gamma$ and $\mu$ ) and the initial condition of both susceptible and infected individuals. Furthermore, it is demonstrated that considering parameter-dependent sets of initial conditions is crucial in some cases.

For the SIR model (which can be regarded as a limit case of the SIRS model with $\mu=0$ ), however, we found that it is observable and identifiable, but not jointly observable-identifiable. Specifically, if both the initial condition and the parameters are unknown, we can only determine $\gamma$ , $\beta/k$ , $\beta S_{0}$ and $kI_{0}$ . This highlights the critical role of the SIRS model’s cyclic dynamics (loss of immunity) in achieving joined observability-identifiability. Based on these findings, we proposed methods to distinguish data from SIR and SIRS models theoretically.

In Section 5, we presented additional epidemiological examples, using one and two observations, for which our framework can also be applied. These examples illustrate the adaptability of the proposed methodologies across different model structures and observational scenarios.

Finally, in Section 6 we demonstrated the applicability of our theoretical results through numerical experiments, using an SIRS model and an SIR model, each under two different scenarios:

•

Scenario 1: The observations are continuous, for which we considered the numerical solutions at each time-step and assumed perfect knowledge of the derivatives of the observations. We solved this scenario through the procedure of Algorithm 4. All tests approximated successfully the exact solution of the SIRS case.
•

Scenario 2: The observations in this scenario were considered to be daily reported (i.e., discrete), and we used Ordinary Least Squares for the calibration of the parameters and presented the results of 5 representative tests for the SIR case. The results indicated that the method reliably captured most of the times the key identifiable combinations ( $\gamma$ , $\beta/k$ , $\beta S(0)$ ) and $\mu$ approximately 0, as expected from the theory exposed in Section 4.2. However, two tests produced approximately the same parameter vectors corresponding to an SIRS model. This result illustrates the difficulty of distinguishing between the two models during early epidemic stages when only a portion of infected individuals is observed, as discussed in Section 4.2.1.

This study demonstrates the feasibility of the proposed methods for parameter estimation and model discrimination in epidemiological systems. Future work should explore extensions to scenarios with more realistic features, such as noisy data (see, e.g., [22]) or bounded measurement ranges (e.g., related to interval observers, see [57]).

Acknowledgments

This work was carried out with financial support from the projects PID2019-106337GB-I00 and PID2023-146754NB-I00, funded by MCIU/AEI/10.13039/501100011033 and FEDER, EU, the project PCI2024-153478, funded by the European M-ERA.Net program, and the French National Research Agency through the ANR project NOCIME (New Observation and Control Issues Motivated by Epidemiology) for the period 2024-2027. In addition, author A.B. Kubik was supported by an FPU predoctoral fellowship and a mobility grant, both provided by the Ministry of Science, Innovation and Universities of the Spanish Government.

References

[1] F. Anstett-Collin, L. Denis-Vidal, and G. Millérioux. A priori identifiability: An overview on definitions and approaches. Annual Reviews in Control, 50:139–149, 2020.
[2] J. Arino and S. Portet. A simple model for COVID-19. Infectious Disease Modelling, 5:309–315, 2020.
[3] G. Bastin and D. Dochain. On-line Estimation and Adaptive Control of Bioreactors. Process Measurement and Control. Elsevier, 1990.
[4] A. Bicchi, D. Prattichizzo, A. Marigo, and A. Balestrino. On the observability of mobile vehicles localization, pages 142–147. World Scientific Pub Co Inc, 1999.
[5] M. Bocher. The theory of linear dependence. Annals of Mathematics, 2(1/4):81–96, 1900.
[6] M. Bongarti, L. Galvan, L. Hatcher, M. Lindstrom, C. Parkinson, C. Wang, and A. Bertozzi. Alternative SIAR models for infectious diseases and applications in the study of non-compliance. Mathematical Models and Methods in Applied Sciences, 32(10):1987–2015, 2022.
[7] F. Brauer, C. Castillo-Chavez, and Z. Feng. Models for Influenza. Mathematical Models in Epidemiology, 69:311–350, 2019.
[8] F. Brauer, C. Castillo-Chavez, A. Mubayi, and S. Towers. Some models for epidemics of vector-transmitted diseases. Infectious Disease Modelling, 1(1):79–87, 2016.
[9] F. Brauer and C. Castillo-Chávez. Mathematical Models in Population Biology and Epidemiology. Springer, 2012.
[10] M. Capistrán, M. Moreles, and B. Lara. Parameter estimation of some epidemic models. The case of recurrent epidemics caused by respiratory syncytial virus. Bulletin of Mathematical Biology, 71(8):1890–1901, 2009.
[11] Y. Chen and J. Zhang. Dynamic analysis and optimal control of influenza model with discrete age structure. Mathematical Methods in the Applied Sciences, 47(6):4260–4282, 2024.
[12] G. Chowell. Fitting dynamic models to epidemic outbreaks with quantified uncertainty: A primer for parameter uncertainty, identifiability, and forecasts. Infectious Disease Modelling, 2(3):379–398, 2017.
[13] N. Cunniffe, F. Hamelin, A. Iggidr, A. Rapaport, and G. Sallet. Observability, Identifiability and Epidemiology – A primer. Springer, 2024.
[14] V. H. D.C. Huong and H. Trinh. Design of event-triggered interval functional observers for systems with input and output disturbances. Mathematical Methods in the Applied Sciences, 44(18):13968–13978, 2021.
[15] L. Denis-Vidal, G. Joly-Blanchard, and C. Noiret. Some effective approaches to check the identifiability of uncontrolled nonlinear systems. Mathematics and Computers in Simulation, 57:35–44, 2001.
[16] S. Diop and M. Fliess. Nonlinear observability, identifiability, and persistent trajectories. In Proceeding of the 30th on Decision and Control, 1981.
[17] A. Dénes and A. Gumel. Modeling the impact of quarantine during an outbreak of Ebola virus disease. Infectious Disease Modelling, 4:12–27, 2019.
[18] N. Evans, L. White, M. Chapman, K. Godfrey, and M. Chappell. The structural identifiability of the susceptible infected recovered model with seasonal forcing. Mathematical Biosciences, 194:175–197, 2005.
[19] M. Fang and P. A. Bliman. Modelling, Analysis, Observability and Identifiability of Epidemic Dynamics with Reinfections. 2022 European Control Conference (ECC), pages 278–283, 2022.
[20] Z. Feng, C. Castillo-Chavez, and A. Capurro. A Model for Tuberculosis with Exogenous Reinfection. Theoretical Population Biology, 57(3):235–247, 2000.
[21] Z. Feng and J. Velasco-Hernández. Competitive exclusion in a vector-host model for the dengue fever. Journal of Mathematical Biology, 35:523–544, 1997.
[22] M. Ferrández, B. Ivorra, J. Redondo, A. Ramos, and P. Ortigosa. A multi-objective approach to identify parameters of compartmental epidemiological models – Application to Ebola Virus Disease epidemics. Communications in Nonlinear Science and Numerical Simulation, 120:107165, 2023.
[23] R. Gani and S. Leach. Transmission potential of smallpox in contemporary populations. Nature, 414:748–751, 2001.
[24] J. Gauthier, H. Hammouri, and S. Othman. A Simple Observer for Nonlinear Systems – Applications to Bioreactors. IEEE Transactions on Automatic Control, 37(6):875–880, 1992.
[25] D. Geller and I. Klein. Angles-Only Navigation State Observability During Orbital Proximity Operations. Journal of Guidance, Control, and Dynamics, 37(6):1976–1983, 2014.
[26] M. Gevers, A. Bazanella, D. Coutinho, and S. Dasgupta. Identifiability and excitation of linearly parametrized rational systems. Automatica, 63:38–46, 2016.
[27] S. Glad and L. Ljung. Parametrization of nonlinear model structures as linear regressions. IFAC Proceedings Volumes, 23(8/3):317–321, 1990.
[28] S. Glad and L. Ljung. On global identifiability for arbitrary model parametrizations. Automatica, 30(2):265–276, 1994.
[29] G. Goodwin and R. Payne. Dynamic System Identification: Experiment Design and Data Analysis. Academic Press, 1977.
[30] F. Hamelin, A. Iggidr, A. Rapaport, G. Sallet, and M. Souza. About the identifiability and observability of the SIR epidemic model with quarantine. IFAC–PapersOnLine, 56(2):4025–4030, 2023.
[31] R. Hermann and A. Krener. Nonlinear Controllability and Observability. IEEE Transactions on Automatic Control, 22(5):728–740, 1977.
[32] S. Holmsen, S. Eidnes, and S. Riemer-Sørensen. Pseudo-hamiltonian system identification. Journal of Computational Dynamics, 11(1):59–91, 2024.
[33] B. Ivorra, M. Ferrández, M. Vela-Pérez, and A. Ramos. Mathematical modeling of the spread of the coronavirus disease 2019 (COVID-19) taking into account the undetected infections. The case of China. Communications in Nonlinear Science and Numerical Simulations, 88:105303, 2020.
[34] R. Jain, S. Narasimhan, and N. Bhatt. A priori parameter identifiability in models with non-rational functions. Automatica, 109:108513, 2019.
[35] G. Jerónimo, M. P. Millán, and P. Solernó. Identifiability from a Few Species for a Class of Biochemical Reaction Networks. Bulletin of Mathematical Biology, 81:2133––2175, 2019.
[36] Q. Jiang, Z. Liu, L. Wang, and R. Tan. A tuberculosis model with early and late latency, imperfect vaccination, and relapse: An application to China. Mathematical Methods in the Applied Sciences, 46(9):10929–10946, 2023.
[37] R. Kalman. A New Approach to Linear Filtering and Prediction Problems. Journal of Basic Engineering, 82(1):35–45, 1960.
[38] R. Kalman. On the General Theory of Control Systems. IFAC Proceedings Volumes, 1(1):491–502, 1960.
[39] S. Kepley and T. Zhang. A constructive proof of the Cauchy–Kovalevskaya theorem for ordinary differential equations. Journal of Fixed Point Theory and Applications, 23(7):1–23, 2021.
[40] W. Kermack and A. McKendrick. A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 115(772):700–721, 1927.
[41] H. Khalil. Nonlinear Systems (3rd Edition). Prentice Hall, 1996.
[42] M. Komatsu and T. Yaguchi. Method for estimating hidden structures determined by unidentifiable state-space models and time-series data based on the gröbner basis. Preprint in arXiv: arXiv.2012.11906, 2020.
[43] M. Kretzschmar, P. Teunis, and R. Pebody. Incidence and reproduction numbers of pertussis: Estimates from serological and social contact data in five European countries. PLoS Medicine, 7(6):e1000291, 2010.
[44] A. Kubik, A. Rapaport, B. Ivorra, and A. Ramos. About identifiability and observability for a class of dynamical systems. Preprint in arXiv: https://doi.org/10.48550/arXiv.2408.11482, 2024.
[45] Y. Liu and J. Rocklöv. The reproductive number of the delta variant of SARS-CoV-2 is far higher compared to the ancestral SARS-CoV-2 virus. Journal of Travel Medicine, 28(7):taab124, 2021.
[46] Y. Liu and J. Rocklöv. The effective reproductive number of the Omicron variant of SARS-CoV-2 is several times relative to Delta. Journal of Travel Medicine, 29(3):taac037, 2022.
[47] D. Luenberger. Observers for Multivariable Systems. IEEE Transactions on Automatic Control, 11(2):190–197, 1966.
[48] A. Martinelli. Observability: A New Theory Based on the Group of Invariance. Society for Industrial and Applied Mathematics, 2020.
[49] G. Massonis, J. Banga, and A. Villaverde. Structural identifiability and observability of compartmental models of the COVID-19 pandemic. Annual Reviews in Control, 51:441–459, 2021.
[50] MathWorks. Help Center – lsqcurvefit. https://www.mathworks.com/help/optim/ug/lsqcurvefit.html.
[51] M. Medvedeva, T. Simos, C. Tsitouras, and V. Katsikis. Direct estimation of SIR model parameters through second-order finite differences. Mathematical Methods in the Applied Sciences, 44(5):3819–3826, 2021.
[52] A. Ovchinnikov, A. Pillay, G. Pogudin, and T. Scanlon. Computing all identifiable functions of parameters for ODE models. Systems & Control Letters, 157:105030, 2021.
[53] A. Ovchinnikov, G. Pogudin, and P. Thompson. Input-output equations and identifiability of linear ODE models. IEEE Transactions on Automatic Control, 68(2):812–824, 2023.
[54] N. Parolini, L. Dede’, G. Ardenghi, and A. Quarteroni. Modelling the COVID-19 epidemic and the vaccination campaign in Italy by the SUIHTER model. Infectious Disease Modelling, 7(2):45–63, 2022.
[55] A. Ramos, M. Ferrández, M. Vela-Pérez, A. Kubik, and B. Ivorra. A simple but complex enough $\theta$ -SIR type model to be used with COVID-19 real data. Application to the case of Italy. Physica D: Nonlinear Phenomena, 421:132839, 2021.
[56] A. Ramos, M. Vela-Pérez, M. Ferrández, A. Kubik, and B. Ivorra. Modeling the impact of SARS-CoV-2 variants and vaccines on the spread of COVID-19. Communications in Nonlinear Science and Numerical Simulation, 102:105937, 2021.
[57] A. Rapaport and D. Dochain. Interval observers for biochemical processes with uncertain kinetics and inputs. Mathematical Biosciences, 193(2):235–253, 2005.
[58] W. Rudin. Principle of Mathematical Analysis (3rd Edition). McGraw Hill, 1976.
[59] M. Saccomani, S. Audoly, and L. D’Angió. Parameter identifiability of nonlinear systems: the role of initial conditions. Automatica, 39:619–632, 2003.
[60] E. Scharbarg, C. H. Moog, N. Mauduit, and C. Califano. From the hospital scale to nationwide: observability and identification of models for the COVID-19 epidemic waves. Annual Reviews in Control, 50:409–416, 2020.
[61] R. Seck, D. Ngom, B. Ivorra, and A. Ramos. An optimal control model to design strategies for reducing the spread of the Ebola virus disease. Mathematical Biosciences and Engineering, 19(2):1746–1774, 2022.
[62] R. Seck, D. Ngom, B. Ivorra, and A. Ramos. An age structured mathematical model for studying Malaria transmission dynamics. Applications to some areas of Senegal. Mathematics and Computers in Simulation, 229:392–408, 2025.
[63] E. Tunali and T.-J. Tarn. New results for identifiability of nonlinear systems. IEEE Transactions on Automatic Control, 32(2):146–154, 1987.
[64] E. Walter and L. Pronzato. Identification of parametric models from experimental data. Springer, 1997.
[65] L. Wang, R. Ortega, and A. Bobtsov. Observability is sufficient for the design of globally exponentially stable state observers for state-affine nonlinear systems. Automatica, 149:110838, 2023.
[66] Z. Wong, C. Bui, A. Chughtai, and C. Macintyre. A systematic review of early modelling studies of Ebola virus disease in West Africa. Epidemiology & Infection, 145(6):1069–1094, 2017.
[67] World Health Organization. Consensus document on the epidemiology of severe acute respiratory syndrome (SARS). https://iris.who.int/handle/10665/70863.
[68] World Health Organization. Ebola disease caused by Sudan ebolavirus – Uganda. https://www.who.int/emergencies/disease-outbreak-news/item/2022-DON423.
[69] World Health Organization. WHO Coronavirus (COVID-19) Dashboard. https://covid19.who.int/.