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

Alicja B Kubik, Benjamin Ivorra, Alain Rapaport, Ángel M Ramos

► To cite this version:

Alicja B Kubik, Benjamin Ivorra, Alain Rapaport, Ángel M Ramos. Identifiability and Observability Analysis for Epidemiological Models: Insights on the SIRS Model. Mathematical Methods in the Applied Sciences, In press, 40 p. 10.1002/mma.70444. hal-05416258

HAL Id: hal-05416258

https://hal.inrae.fr/hal-05416258v1

Submitted on 15 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Identifiability and Observability Analysis for Epidemiological
Models: Insights on the SIRS Model

Alicja B. Kubik*¹, Benjamin Ivorra†², Alain Rapaport‡³, and Ángel M. Ramos§²

¹National Laboratory for Health Security, Bolyai Institute, Univ. of Szeged, 6720 Szeged, Hungary

²MOMAT, Instituto de Matemática Interdisciplinar, Univ. Complutense de Madrid, 28040 Madrid, Spain

³MISTEA, Univ. Montpellier, INRAE, Institut Agro, 34060 Montpellier, France

Abstract

The problems of observability and identifiability have been of great interest as previous steps to estimat-
ing 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 non-
linear systems, including non-analytic systems have received far less attention and, to the best of our
knowledge, no unified constructive methodology has been proposed to assess and guarantee parameter
and state recoverability in this setting. Some symbolic tools provide automated analyses for rational or
nonlinear, analytical systems, offering qualitative identifiability and observability verdicts that are some-
times incomplete. In this work, we introduce a family of efficient and fully constructive procedures that
can enable explicit recovery of the unknown parameters and/or initial conditions, whenever possible, for
a large class, not necessarily rational or analytic, nonlinear ODE systems. 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
joint 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; Joint observability-identifiability; Parameter identifiability; State ob-servability; 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

<footnote>*akubik@ucm.es</footnote>

<footnote>†ivorra@mat.ucm.es</footnote>

<footnote>‡alain.rapaport@inrae.fr</footnote>

<footnote>§angel@mat.ucm.es</footnote>
up in the early 1900s (see [48]), and have been regularly used since then to estimate the dynamics of different diseases (for example, influenza [9, 15], tuberculosis [28, 44], Ebola [22, 73], or vector-transmitted diseases such as malaria [72], dengue fever or the Zika virus [11, 29]), with a significant surge in applications in recent years due to the COVID-19 pandemic (see, for example, [2, 8, 41, 64, 67]). 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, comparatively much less 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 instance, [13, 16, 37, 20, 71, 61, 76]). 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* [46]). 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 uniquely 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 [19]) 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 [39]); 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. [58]). 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 [56], high-gains observers [32], or
the Kalman filter [45]. 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 [19, 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 [4, 6, 33, 40, 42, 77]). In particular, the literature about joint observability-identifiability of epidemiological models has increased in recent years (see [19, 27, 55, 38, 59] and the references therein).

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., [23, 70]). While there are fields where one can perform
different experiments with chosen initial conditions in order to identify the parameters (see, e.g., [23, 43]), 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 [52], 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., [5, 24, 25, 34, 35, 36]), and the study of
nonlinear, non-rational systems is scarce (see [1, 3]). For rational systems, a typical goal is finding the input-
output (IO) equations in order to determine if the parameters are identifiable. These are structural equations
of the system under study computed through a (classical) “elimination” step which consists in finding an
expression that relates the outputs and their derivatives, the inputs and their derivatives (if any), and the
parameters, independently of the state variables. However, we consider that in these works there lacks a
direct, unified argument to conclude that the parameters can be effectively recovered from these equations.
In [25], the authors perform a second step which aims to answer this question. The authors do it using
probabilistic search, which turns out to be often efficient but that could miss some isolated singular points.
For analytic nonlinear systems, we can find the software STRIKE-GOLDD developed in [75], which assesses
on local identifiability relying on a generalization of the ORC aforementioned, typically requiring high-order
derivatives. In [52], general nonlinearities, not necessarily analytic, 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 joint observability-identifiability in nonlinear systems
of ODEs.

In [52], we addressed several of the mentioned gaps; in particular, the lacking argument ensuring recover-
ability of unknowns 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 [23] 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 [62] 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 [52], 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 suited 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 joint 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 \ge 0$, represented by $y$:

$$ \begin{cases} \dot{S} = -\beta SI + \mu R, \\ \dot{I} = \beta SI - \gamma I, \\ \dot{R} = \gamma I - \mu R, \\ y = kI, \end{cases} \quad \forall t \ge 0. \quad (1) $$

Here, $\beta$, $\gamma$, and $\mu$ denote respectively the transmission, recovery, and immunity-loss rates. This model is particularly relevant for understanding the dynamics of diseases such as influenza or 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 instance, [66]).

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 [52]. The presented SIRS model with output $y = kI$ fits within the class of systems described in Section 3 and serves as a motivating example. Indeed, in Section 4 we make use of this approach to prove its joint 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$, under 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 joint observability-identifiability, highlighting a key distinction between these two models that allows us to distinguish them at 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 [52] for analyzing observability, identifiability and joint 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 joint observability-identifiability. Unlike in [52], 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 [52], 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 parameters 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

$$ egin{cases} \frac{dx}{dt}(t; \xi, \theta) = f(x(t; \xi, \theta), \theta), & x(0; \xi, \theta) = \xi, \\ y_{(\xi, \theta)}(t) = h(x(t; \xi, \theta), \theta), \end{cases} \qquad (2) $$

where $f(\cdot, \cdot) : \Omega \times \Theta \to \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} \to \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 \to \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 $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 \in \Theta$ produces a different output at some time $t \in S$, i.e.,

$$h(x(t; \xi, \theta_0), \theta_0) \neq h(x(t; \xi, \theta'_0), \theta'_0).$$

Equivalently, if $h(x(t; \xi, \theta_0), \theta_0) = h(x(t; \xi, \theta'_0), \theta'_0)$, for all $t \in S$ and any $\xi \in \Omega_{\theta_0}$, then $\theta_0 = \theta'_0$. Moreover, if $\Omega_{\theta} = \Omega$, for each $\theta \in \Theta$, we say that System (2) is identifiable on $\Theta$ in $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 parameters 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 $S \subset \mathcal{I}$ with parameters in $\Theta$ whether, for any $\theta \in \Theta$, any two different $x_0, x'_0 \in \Omega$ produce a different output at some time $t \in S$, i.e.,

$$h(x(t; x_0, \theta), \theta) \neq h(x(t; x'_0, \theta), \theta).$$

Equivalently, if $h(x(t; x_0, \theta), \theta) = h(x(t; x'_0, \theta), \theta)$, for all $t \in S$ and any $\theta \in \Theta$, then $x_0 = x'_0$.

**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 $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 $S = \mathcal{I}$.

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

**Example 1.** Consider the following system:

$$\begin{cases} x'(t) = k, & x(0) = \xi, \\ y_{(\xi,k)}(t) = \max(1, x(t)), \end{cases} \quad (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 \ge 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 \ge 1/\min\{k_1, k_2\}$. However, if $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 S$, and hence the system is not identifiable on $\Theta$ in $S$ with initial conditions in $\Omega$. $\square$

**Example 2.** We consider the same case as the one shown in Example 1. 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 \in \Omega$, we have that

$$\max(1, kt + x_0) = kt + x_0 \neq kt + x'_0 = \max(1, kt + x'_0),$$
for all $t \ge 1/k$. However, if $S = [0, 1/2] \subset \mathcal{I}$ and $x_0, x'_0 \in [0, 1/2] \subset \Omega$, with $x_0 \ne x'_0$, then

$$ \max(1, kt + x_0) = 1 = \max(1, kt + x'_0), $$

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

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 joint observability-identifiability of a system (see [19, 74]).

We are now going to consider that we need to recover both $x_0$ and $\theta_0$.

**Definition 5 (Joint observability-identifiability in a time set).** System (2) is jointly observable-identifiable on $\Gamma_{\Theta}$ in $S \subset \mathcal{I}$ if, for any $(x_0, \theta_0) \in \Gamma_{\Theta}$, a different $(x'_0, \theta'_0) \in \Omega \times \Theta$ produces a different output at some time $t \in S$, i.e.,

$$ h(x(t; x_0, \theta_0), \theta_0) \neq h(x(t; x'_0, \theta'_0), \theta'_0). $$

Equivalently, if $h(x(t; x_0, \theta_0), \theta_0) = h(x(t; x'_0, \theta'_0), \theta'_0)$ for all $t \in S$, this implies that $x_0 = x'_0$ and $\theta_0 = \theta'_0$.

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

Furthermore, joint observability-identifiability in a time set $S$ implies joint observability-identifiability in any time set $\mathcal{A}$ such that $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 joint 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 joint 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

For completeness, we briefly recall the well-known concept of the Lie derivative and the way it is computed, since it will be used repeatedly throughout the paper. This concept is fundamental for identifiability and observability analysis, as it provides a way to express the time derivatives of the output as functions of the state variables and parameters of the system.

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

$$ \mathcal{L}_{F, \varphi, d}: \begin{array}{rcl} \Omega & \to & \mathbb{R}^{m \times (d+1)} \\ \xi & \mapsto & (L_F^0 \varphi(\xi), L_F^1 \varphi(\xi), \dots, L_F^d \varphi(\xi)), \end{array} $$

where $L_F^0 \varphi = \varphi$ and, if $d \ge 1$, $L_F^1 \varphi \in 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) = D\varphi(\xi)F(\xi), $$
where Dφ is the $m \times n$ Jacobian matrix of φ. In particular, if $z(\cdot; \xi)$ is a solution of

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

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

$$ L_F^1 \varphi(z(t; \xi)) = D\varphi(z(t; \xi))F(z(t; \xi)) = \frac{d}{dt}\varphi(z(t; \xi)). \qquad (5) $$

Then,

$$ L_F^1 \varphi(\xi) = \left. \frac{d}{dt} \varphi(z(t; \xi)) \right|_{t=0}. \qquad (6) $$

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

$$ L_F^k \varphi = L_F^1 (L_F^{k-1} \varphi). $$

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

$$ L_F^k \varphi(z(t; \xi)) = \left. \frac{d^k}{dt^k} \varphi(z(t; \xi)) \right|_{t=0} \quad \text{and, in particular, } \quad L_F^k \varphi(\xi) = \left. \frac{d^k}{dt^k} \varphi(z(t; \xi)) \right|_{t=0}. $$

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

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

Then, for $k$, we have

$$ \frac{d^k}{dt^k} \varphi(z(t; \xi)) = \left. \frac{d^k}{dt^k} L_F^{k-1} \varphi(z(t; \xi)) \right|_{t=0} = D L_F^{k-1} \varphi(z(t; \xi)) F(z(t; \xi)) = L_F^k \varphi(z(t; \xi)), \quad \text{and} \quad L_F^k \varphi(\xi) = \left. \frac{d^k}{dt^k} \varphi(z(t; \xi)) \right|_{t=0}, \\ \text{as we wanted to prove.} \qquad \square $$

Let us now compute $\mathcal{L}_{F,\phi,d}$ for System (2). Given $\theta \in \Theta$, denote $h_\theta(\cdot) = h(\cdot, \theta)$ and $f_\theta(\cdot) = f(\cdot, \theta)$, and assume $h_\theta \in C^d(\Omega; \mathbb{R}^m)$, $f_\theta \in C^{\max\{0,d-1\}}(\Omega; \mathbb{R}^n)$, for some $d \in \mathbb{N} \cup \{0\}$. Then, $y_{(\xi,\theta)} \in C^d(\mathcal{I}; \mathbb{R}^m)$, since

$$ \dot{y}_{(\xi,\theta)}(t) = \left. \frac{d}{dt} h(x(t; \xi, \theta)) \right|_{t=0} = Dh_{\theta}(x(t; \xi, \theta)) f_{\theta}(x(t; \xi, \theta)) \implies \dot{y}_{(\xi,\theta)}(t) = L_{f_{\theta}}^{1} h_{\theta}(x(t; \xi, \theta)), \quad \forall t \ge 0. $$

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

$$ y_{(\xi,\theta)}^{(k)}(t) = \left. \frac{d^k}{dt^k} h(x(t; \xi, \theta)) \right|_{t=0} = L_{f_{\theta}}^k h_{\theta}(x(t; \xi, \theta)) \quad \text{and, hence, } y_{(\xi,\theta)}^{(k)}(0) = \left. \frac{d^k}{dt^k} h(x(t; \xi, \theta)) \right|_{t=0} = L_{f_{\theta}}^k h_{\theta}(\xi), $$

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

$$ \mathcal{L}_{f_\theta, h_\theta, d}(\xi) = (y_{(\xi,\theta)}^{(0)}(0), \dots, y_{(\xi,\theta)}^{(d)}(0)). $$

This property underlies most observability and identifiability analyses, as it allows reconstructing states or parameters from instantaneous measurable outputs and their derivatives.

In the following, we denote $y_{(\xi,\theta)}$, $\dot{y}_{(\xi,\theta)}$ and $y_{(\xi,\theta)}^{(k)}$ as $y$, $\dot{y}$ and $y^{(k)}$, respectively, when the context is clear. We also let $h_\theta = (h_{\theta,1}, \dots, h_{\theta,m})$, where $m \in \mathbb{N}$ is the number of scalar outputs.

We now proceed to present the main theoretical results and corresponding algorithms for recovering both parameters and initial conditions.
## 3.3 Main results

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

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

**Theorem 1.** Let $\Omega \subset \mathbb{R}^n$ be positively invariant with respect to the ODE system given in (2), $h_{\theta,i} \in C^{d_i}(\Omega; \mathbb{R})$ for some $d_i \in \mathbb{N} \cup \{0\}$, $i \in \{1, \dots, m\}$, and $f_{\theta} \in C^{d-1}(\Omega; \mathbb{R}^n)$, where $d = \max\{1, d_1, \dots, d_m\}$, for any $\theta \in \Theta$. If

$$ \mathcal{L}_{f_{\theta},h_{\theta},\{d_1,\dots,d_m\}} : \begin{array}{rcl} \Omega & \rightarrow & \mathbb{R}^{m+d_1+\dots+d_m} \\ \xi & \mapsto & (\mathcal{L}_{f_{\theta},h_{\theta,1},d_1}(\xi),\dots,\mathcal{L}_{f_{\theta},h_{\theta,m},d_m}(\xi)) \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 \in \Omega$ such that

$$ h(x(t; x_0, \theta), \theta) = h(x(t; x'_0, \theta), \theta), \quad \forall t \in [a, b). $$

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

$$ h(x(t; \xi, \theta), \theta) = h(x(t; \xi', \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_{(\xi,\theta),i}^{(k)}(t) = y_{(\xi',\theta),i}^{(k)}(t), \quad \forall t \in [0, b - \tilde{t}), \quad k \in \{0, \dots, d_i\}, \quad i \in \{1, \dots, m\}, $$

where $y_{(\xi,\theta)}^{(k)} = (y_{(\xi,\theta),1}^{(k)}, \dots, y_{(\xi,\theta),m}^{(k)})$. In particular, $y_{(\xi,\theta),i}^{(k)}(0) = y_{(\xi',\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'). $$

Since $\mathcal{L}_{f_{\theta},h_{\theta},\{d_1,\dots,d_m\}}$ is injective in $\Omega$, for any $\theta \in \Theta$, and $\xi, \xi' \in \Omega$ due to the positive invariance, this implies that $\xi = \xi'$, i.e., $x(\tilde{t}; x_0, \theta) = x(\tilde{t}; x'_0, \theta)$. Due to the uniqueness of solutions of System (2), this implies that $x_0 = x'_0$.

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

Theorem 1 states that a system is observable if the mapping $\mathcal{L}_{f_{\theta},h_{\theta},\{d_1,\dots,d_m\}}$, which collects the outputs and their successive Lie derivatives up to some finite order, is injective in the state space (or some positively invariant subset). Intuitively, each derivative of the measured output provides new information about the hidden variables; if no two distinct states can produce the same sequence of derivatives, the state of the system can be uniquely reconstructed from the observations. This can be viewed as the nonlinear analogue of the observability rank condition in linear observability theory.

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) := (y_{(\xi,\theta),1}^{(0)}(t), \dots, y_{(\xi,\theta),1}^{(d_1)}(t), \dots, y_{(\xi,\theta),m}^{(0)}(t), \dots, y_{(\xi,\theta),m}^{(d_m)}(t)), $$
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 [52], which we now adapt to Definition 1 and provides a sufficient condition to ensure recoverability of parameters.

**Theorem 2.** Let $h_{\theta,i} \in \mathcal{C}^{d'_i}(\Omega; \mathbb{R})$, for some $d'_i \in \mathbb{N}_0^1$, $i \in \{1, \dots, m\}$, 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$. Consider $\mathcal{D} \subset \mathbb{R}^{d'_1+\dots+d'_{m}+m}$ such that, for all $(\xi, \theta) \in \Gamma_{\Theta}$, $\mathbf{y}^{(d_1, \dots, d_m)}(t) \in \mathcal{D}$, for all $t \in \mathcal{I}$. Let $\mathcal{S} \subset \mathcal{I}$ be such that every connected component of $\mathcal{S}$ contains an open interval. Assume there exist $g: \mathcal{D} \to \mathbb{R}^{q+p}$ and $r: \Theta \to \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, \dots, d'_m)}) = \sum_{l=1}^{q_j} r_{j,l}(\theta) g_{j,l}(\mathbf{y}^{(d'_1, \dots, d'_m)}), \quad (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, \dots, d'_m)}(t))\}_{l=1}^{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 \in \Theta$ such that

$$h_{\theta_0}(x(t; \xi, \theta_0)) = h_{\theta'_0}(x(t; \xi, \theta'_0)), \quad \forall t \in \mathcal{S},$$

i.e.

$$y_{(\xi, \theta_0)}(t) = y_{(\xi, \theta'_0)}(t), \quad \forall t \in \mathcal{S}.$$

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

$$y_{(\xi, \theta_0),i}^{(k)}(t) = y_{(\xi, \theta'_0),i}^{(k)}(t), \quad \forall t \in \mathcal{S}, k \in \{0, \dots, d'_i\}, i \in \{1, \dots, m\}.$$

Since

$$g_{j,0}(\mathbf{y}^{(d'_1, \dots, d'_m)}_{(\xi, \theta_0)}) - g_{j,0}(\mathbf{y}^{(d'_1, \dots, d'_m)}_{(\xi, \theta'_0)}) \equiv 0$$

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

$$\sum_{l=1}^{q_j} (r_{j,l}(\theta_0) - r_{j,l}(\theta'_0)) g_{j,l}(\mathbf{y}^{(d'_1, \dots, d'_m)}_{(\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, \dots, d'_m)}_{(\xi, \theta_0)})\}_{l=1}^{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)$. Since $r$ is injective, we have $\theta_0 = \theta'_0$. 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}$. $\square$

Theorem 2 formalizes how the model outputs and their derivatives can determine identifiability of the system. The key idea is to separate the time dependence, represented by the functions $g_{j,l}$, from the parameter dependence, encoded by the coefficients $r_{j,l}(\theta)$. If $\{g_{j,l}(\mathbf{y}^{(d'_1, \dots, d'_m)}_{(\xi, \theta_0)})\}_{l=1}^{q_j}$, for any $j \in \{1, \dots, p\}$, are linearly independent over time and the map $r$ is injective, then we can ensure that all outputs are generated by a unique parameters vector.

<footnote n="1">¹We use a different notation with respect to Theorem 1 to explicitly state that these orders may be different.</footnote>
**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 1, the conditions presented in Theorems 1 and 2 will allow us to have joint observability-identifiability of System (2). In [23, 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.

**Remark 2.** For rational systems, the equations obtained from System (7) are equivalent to the so-called input-output equations produced by the DAISY algorithm [5], which derives differential polynomial equations linking the model outputs and their derivatives, the model inputs (if any) and their derivatives, and the unknown parameters using characteristic sets and Gröbner bases. Our formulation is equivalent modulus a convenient rewriting that allows us to verify conditions (C2) and (C3). We will illustrate the importance of this rewriting throughout Sections 4 and 5. On the other hand, the software StructuralIdentifiability.jl [25] concludes the recoverability of parameters using a probabilistic search. However, this search may miss particular states and parameters that give linear dependence, which will be as well important in the previously mentioned sections.

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 joint observability-identifiability of our model and recover the initial condition and parameters 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 $S \subset I$, $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 $S$ and we can reconstruct the pair $(x_0, \theta_0)$ uniquely using the values of $y_{(x_0, \theta_0)}$ and its derivatives at, at most, $q+1 = q_1 + \cdots + q_p + 1$ suitable values of $t \in S$.

*Proof.* Let $\phi_{j,l}(t) = g_{j,l}(\mathbf{y}_{(x_0, \theta_0)}^{(d'_1, \dots, d'_m)}(t))$, $t \in 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 $S$, for each $j \in \{1, \dots, p\}$. Then, for every $j \in \{1, \dots, p\}$, there exist (see [52, Lemma 1]) $q_j$ different time instants $\{t_{j,\ell}\}_{\ell=1}^{q_j} \subset 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

$$ (\phi_{j,1}(t_{j,\ell}) \cdots \phi_{j,q_j}(t_{j,\ell})) \sigma = \phi_{j,0}(t_{j,\ell}), \quad \forall \ell \in \{1, \dots, q_j\}. \quad (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) \to \Theta$, we can recover our original parameters vector

$$ \theta_0 = r^{-1}(\sigma_{1,1}, \dots, \sigma_{p,q_p}). $$

Finally, to recover $x_0$, take some time $\tilde{t} \in 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}_{f_{\theta_0}, h_{\theta_0}, \{d_1, \dots, d_m\}}^{-1}(\mathbf{y}_{(x_0, \theta_0)}^{(d_1, \dots, d_m)}(\tilde{t})), $$

noticing that $\mathbf{y}_{(x_0, \theta_0)}^{(d_1, \dots, d_m)}(\tilde{t}) = \mathbf{y}_{(\tilde{\xi}, \theta_0)}^{(d_1, \dots, d_m)}(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 S$, we directly choose $\tilde{t} = 0$.
This is, we have recovered $\theta_0$ and $x_0$ from the data uniquely knowing $y_{(x_0,\theta_0)}$ in $S$, using its value and the value of its derivatives at $q+1$ (at most) different time instants. $\square$

In Algorithm 1 below, we present a procedure to recover the unknown quantities from given observations, following the proof of Theorem 3, and emphasizing the importance of selecting appropriate sets for the initial conditions.

**Algorithm 1.** Assume that we know $y_{(x_0,\theta_0)}$ (satisfying System (2)) in $S \subset \mathcal{I}$, such that every connected of $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}_{(\xi, \theta)}^{(d_1, \dots, d_m)}(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} \to \mathbb{R}^{q+p}$ and $r: \Theta \to \mathbb{R}^p$, for some suitable $\mathcal{D} \subset \mathbb{R}^{d_1'+\dots+d_m'+m}$, $d_1', \dots, d_m' \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 $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 S$ such that

$$ \det(\phi_{j,l}(t_{j,l})_{l,l=1,\dots,q_j}) \neq 0, $$

where $\phi_{j,l} = g_{j,l}(\mathbf{y}_{(x_0,\theta_0)}^{(d_1',\dots,d_m')})$.

Note: These time instants exist due to condition (C3) of Theorem 2, which is verified in Step 2, and [52, 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), \quad (9) $$

where $\phi_{j,0} = g_{j,0}(\mathbf{y}_{(x_0,\theta_0)}^{(d_1',\dots,d_m')})$.

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 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}_{(x_0, \theta_0)}^{(d'_1, \dots, d'_m)}(\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 $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 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 2.

**Remark 3.** Notice that, in order to be able to recover $(x_0, \theta_0)$ following the procedure in Step 4 of Algorithm 1, for each set $\{\phi_{j,l}\}_{l=1}^{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 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 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 1.

**Note:** Although the idea of selecting a suitable number of time instants was already mentioned in [50] and [70, Remark 3], no formal proof nor method to choose these times were provided.

**Lemma 2.** Assume the hypotheses of Theorem 2 are satisfied for some $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, \dots, d'_m)})$ are analytic on $\mathcal{I}$. Then, System (2) is identifiable on $\Theta$ in any semi-open interval $[a, b) \subset S$ with initial conditions consistent with the family $\{\Omega_\theta\}_{\theta \in \Theta}$. If $S$ is connected, it is sufficient for $\{g_{j,l}(\mathbf{y}^{(d'_1, \dots, d'_m)})\}_{\substack{j=1, \dots, p \\ l=1, \dots, q_j}}$ to be analytic on $S$.

*Proof.* Given $\theta_0 \in \Theta$, $\xi \in \Omega_{\theta_0}$, let $\theta'_0 \in \Theta$ such that

$$ h_{\theta_0}(x(t; \xi, \theta_0)) = h_{\theta'_0}(x(t; \xi, \theta'_0)), \quad \forall t \in [a, b), $$

i.e.,

$$ y_{(\xi, \theta_0)}(t) = y_{(\xi, \theta'_0)}(t), \quad \forall t \in [a, b). $$

Then,

$$ y_{(\xi, \theta_0), i}^{(k)}(t) = y_{(\xi, \theta'_0), i}^{(k)}(t), \quad \forall t \in [a, b), $$

$k \in \{0, \dots, d'_i\}$, $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, \dots, d'_m)})$ are analytic in $\mathcal{I}$, then the function

$$ G_{j,\theta}(\mathbf{y}^{(d'_1, \dots, d'_m)}) := \sum_{l=1}^{q_j} r_{j,l}(\theta) g_{j,l}(\mathbf{y}^{(d'_1, \dots, d'_m)}) $$

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

$$ G_{j,\theta}(\mathbf{y}^{(d'_1, \dots, d'_m)}) = g_{j,0}(\mathbf{y}^{(d'_1, \dots, d'_m)}), $$
in $[a, b)$, for any $\theta \in \Theta$. Since

$$g_{j,0}(\mathbf{y}_{(\xi,\theta_0)}^{(d'_1, \dots, d'_m)}) = g_{j,0}(\mathbf{y}_{(\xi,\theta'_0)}^{(d'_1, \dots, d'_m)}),$$

in $[a, b)$, we have

$$G_{j,\theta_0}(\mathbf{y}_{(\xi,\theta_0)}^{(d'_1, \dots, d'_m)}) = G_{j,\theta'_0}(\mathbf{y}_{(\xi,\theta'_0)}^{(d'_1, \dots, d'_m)}),$$

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

$$R_j := \sum_{l=1}^{q_j} (r_{j,l}(\theta_0) - r_{j,l}(\theta'_0)) g_{j,l}(\mathbf{y}_{(\xi,\theta_0)}^{(d'_1, \dots, d'_m)}) \equiv 0$$

in $[a, b)$. Due to the analyticity of $\{g_{j,l}(\mathbf{y}^{(d'_1, \dots, d'_m)})\}_{l=1}^{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}$ ([69, Theorem 8.5]). Therefore, since $\{g_{j,l}(\mathbf{y}^{(d'_1, \dots, d'_m)})\}_{l=1}^{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)$. Since $r$ is injective, this implies $\theta_0 = \theta'_0$.

If $\mathcal{S}$ is connected, $\{g_{j,l}(\mathbf{y}^{(d'_1, \dots, d'_m)})\}_{l=1}^{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, \dots, d'_m)})\}_{l=1}^{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}$. $\square$

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)$ uniquely using the values of $y_{(x_0,\theta_0)}$ and its derivatives at a finite amount of suitable time 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, \dots, d'_m)})\}_{l=1}^{q_j}$ are linearly independent in some $\mathcal{S} \subset \mathcal{I}$ and analytic in $\mathcal{I}$, they are also linearly independent in $[a, b) \subset \mathcal{S}$. Given $j \in \{1, \dots, p\}$, assume that $\{g_{j,l}(\mathbf{y}^{(d'_1, \dots, d'_m)})\}_{l=1}^{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, \dots, d'_m)}(t)) = 0, \forall t \in [a, b).$$

Since $\{g_{j,l}(\mathbf{y}^{(d'_1, \dots, d'_m)})\}_{l=1}^{q_j}$ are analytic in $\mathcal{I}$, then so it is $G_j$. This implies that (see [69, Theorem 8.5]) $G_j \equiv 0$ in $\mathcal{I}$ and, thus, $\{g_{j,l}(\mathbf{y}^{(d'_1, \dots, d'_m)})\}_{l=1}^{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, \dots, d'_m)})\}_{l=1}^{q_j}$ to be 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 [69, Theorem 8.5]) that $G_j \equiv 0$ in $\mathcal{S}$, which leads to the same contradiction. Thus, 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}_{(x_0,\theta_0)}^{(d'_1, \dots, d'_m)}), l \in \{1, \dots, q_j\}$, are linearly independent in $[a, b)$ and we can conclude analogously to Theorem 3, along with Remark 3, 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)$ uniquely 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$ distinct suitable time points. $\square$
We present in Algorithm 2 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 1, 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 1, 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}$. Find maps $g : \mathcal{D} \to \mathbb{R}^{q+p}$ and $r : \Theta \to \mathbb{R}^p$, for some suitable $\mathcal{D} \subset \mathbb{R}^{d'_1+\dots+d'_m+m}$, $d'_1, \dots, d'_m \in \mathbb{N} \cup \{0\}$, $q, p \in \mathbb{N}$, such that (C1), (C2) and (C3) of Theorem 2 are satisfied, being $S \subset \mathcal{I}$ in (C2) and (C3) a connected set satisfying $[a,b) \subset S$, and such that functions $\{g_{j,l}(\mathbf{y}^{(d'_1,\dots,d'_m)})\}_{\substack{j=1,\dots,p \\ l=1,\dots,q_j}}$ are analytic in $S$, for any $(\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 S$ with initial conditions consistent with the family $\{\Omega_{2,\theta}\}_{\theta \in \Theta}$.

Step 4. Same as in Algorithm 1 but such that, for each $j \in \{1, \dots, p\}$, $\{t_{j,l}\}_{l=1}^{q_j} \subset [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 [52, Lemma 1].

Step 7. Same as in Algorithm 1, but such that $\tilde{t} \in \{t_{1,1}, \dots, t_{p,q_p}\}$.

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 4.** Recall that, as in Remark 3, 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 in Algorithm 2, where the needed order of the derivatives is given by the same algorithm.

The use of these algorithms will be illustrated in detail in Section 4.

## 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, [34, 35, 63]). In the methodology presented in Lemma 4 and Algorithm 3 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 I$, and the hypotheses of Theorem 1 and Lemma 2 are satisfied for some $S \subset I$ such that $[a, b) \subset S$. If we further assume that $h_{\theta, i} \in C^{d'_i+\bar{q}-1}(\Omega; \mathbb{R})$, $i \in \{1, \dots, m\}$, $f_{\theta} \in C^{d'+\bar{q}-2}(\Omega; \mathbb{R}^n)$, for any $\theta \in \Theta$, being $\bar{q} = \max\{q_1, \dots, q_p\}$, then System (2) is jointly observable-identifiable and we can reconstruct $(x_0, \theta_0)$ uniquely 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}(Y_{(x_0, \theta_0)}^{(d'_1, \dots, d'_m)})$ 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

$$ \frac{d^k \phi_{j,0}}{dt^k}(t) = \sum_{l=1}^{q_j} r_{j,l}(\theta_0) \frac{d^k \phi_{j,l}}{dt^k}(t), \quad \forall t \in [a, b), \quad k \in \{0, \dots, q_j - 1\}, $$

By the proof of [52, Lemma 3], we know that $\{\phi_{j,l}\}_{l=1}^{q_j}$ are linearly independent in $[a, b)$, and hence there exists some $t_j \in [a, b)$ (see [7]) such that

$$ W_j(t_j) = \det \left( \left( \frac{d^k \phi_{j,l}(t_j)}{dt^k} \right)_{\substack{l=1, \dots, q_j \\ k=0, \dots, q_j-1}} \right) \neq 0, $$

where $W_j$ is the Wronskian of $\{\phi_{j,l}\}_{l=1}^{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) \\ \frac{d\phi_{j,1}}{dt}(t_j) & \cdots & \frac{d\phi_{j,q_j}}{dt}(t_j) \\ \vdots & \ddots & \vdots \\ \frac{d^{q_j-1}\phi_{j,1}}{dt^{q_j-1}}(t_j) & \cdots & \frac{d^{q_j-1}\phi_{j,q_j}}{dt^{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) \\ \frac{d\phi_{j,0}}{dt}(t_j) \\ \vdots \\ \frac{d^{q_j-1}\phi_{j,0}}{dt^{q_j-1}}(t_j) \end{array} \right). $$

Moreover, $W_j$ is also analytic, and therefore it can only vanish at isolated points [69, 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 [52, Theorem 3], we may recover our original parameters 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 [52, Lemma 3].

Hence, we are able to recover $(x_0, \theta_0)$ uniquely 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)$. $\square$

Therefore, if we ask for higher regularity of the components of $y_{(\xi, \theta)}$, $(\xi, \theta) \in \Gamma_{\Theta}$, than required in 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\} \ge 1$ and $q = q_1 + \dots + q_p \ge p$ (see Remark 3). We present next the associated algorithm.

**Algorithm 3.** In a way similar to those of Algorithm 1 and Algorithm 2, 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 1 and Algorithm 2. We only state here the steps that are different from those in Algorithms 1 and 2:

Step 2. Same as in Step 2 of Algorithm 2, along with $y_{(\xi, \theta), i} \in C^{d'_i+\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,l} = g_{j,l}(\mathbf{y}_{(x_0, \theta_0)}^{(d'_1, \dots, d'_m)}), \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( \frac{d^k \phi_{j,l}}{dt^k}(\tilde{t}) \right)_{\substack{l=1, \dots, q_j \\ k=0, \dots, q_j-1}} \right) \neq 0, \quad \forall j \in \{1, \dots, p\}, \qquad (10) $$

where $W_j$ is the Wronskian of $\{\phi_{j,l}\}_{l=1}^{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}) \\ \frac{d\phi_{j,1}}{dt}(\tilde{t}) & \cdots & \frac{d\phi_{j,q_j}}{dt}(\tilde{t}) \\ \vdots & \ddots & \vdots \\ \frac{d^{q_j-1}\phi_{j,1}}{dt^{q_j-1}}(\tilde{t}) & \cdots & \frac{d^{q_j-1}\phi_{j,q_j}}{dt^{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}) \\ \frac{d\phi_{j,0}}{dt}(\tilde{t}) \\ \vdots \\ \frac{d^{q_j-1}\phi_{j,0}}{dt^{q_j-1}}(\tilde{t}) \end{array} \right). \qquad (11) $$

Step 7. Same as in Step 6 of Algorithm 1, but considering that $r(\theta_0)$ is solution of System (11).

Step 8. Same as in Step 7 of Algorithm 2, 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 1.
**Remark 5.** Notice that, as in Remark 3 and Remark 4, 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 + \tilde{q} - 1\} - \text{th derivative}$) at one time such that Step 5 in Algorithm 3 is satisfied. Moreover, if $0 \in 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 3 is satisfied and not need to integrate backwards.

In practice, the recovery processes in the constructive algorithms presented above follow three main steps: (i) determine whether the system is observable by checking the injectivity of a mapping of Lie derivatives; (ii) verify identifiability through the linear relations in (7), ensuring that the time-dependent right-hand side functions are linearly independent; and (iii) use a finite set of outputs and derivative values at selected times to solve explicit linear systems and reconstruct the unknown parameters and/or initial conditions. Algorithms 1–3 correspond to different ways of performing these steps, depending on whether multiple sampling instants or higher-order derivatives are available.

**Remark 6.** Once obtained the parameters vector and initial conditions, 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, \tilde{t}\}_{j=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. [45]). 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 parameters vector. Besides, if we have some analyticity properties, we can use almost any time $t \in 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 [47]), 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 7.** We can apply the proposed algorithms to some systems with piece-wise constant parameters. If the parameters 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 1, 2 and 3 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.

## 3.5 Summary of the proposed methodology

For clarity, we provide here a concise and non-technical summary of the proposed framework to assess joint observability-identifiability in nonlinear dynamical systems.

**Model proposal.** Consider a suitable ODE model along with measurable outputs:

$$ \left\{ \begin{array}{l} \frac{dx}{dt}(t; \xi, \theta) = f(x(t; \xi, \theta), \theta), \quad x(0; \xi, \theta) = \xi, \\ y_{(\xi, \theta)}(t) = h(x(t; \xi, \theta), \theta), \end{array} \right. \qquad (12) $$
where $x$ denotes the state variables, $\theta \in \Theta$ the unknown parameters, and $y$ the measurable outputs. Identify an appropriate initial positively invariant domain $\Omega$ in which the states evolve.

**Find a suitable observability map.** Define a mapping $\mathcal{L}_{f_{\theta},h_{\theta},\{d_1,...,d_m\}}$ that collects $y$ and its derivatives up to some order. If, for any chosen parameters vector in $\Theta$, it is injective in $\Omega$ or a positively invariant subset, observability of the initial conditions is ensured.

**Find structural equations of the system.** Formulate relations of the type

$$g_{j,0}(\mathbf{y}^{(d'_1,...,d'_m)}) = \sum_{l=1}^{q_j} r_{j,l}(\theta)g_{j,l}(\mathbf{y}^{(d'_1,...,d'_m)}),$$

where $g_{j,l}$ are functions of the measured outputs and their derivatives, and $r_{j,l}(\theta)$ depend only on the parameters.

**Identifiability conditions.** If the functions $g_{j,l}(\mathbf{y}^{(d'_1,...,d'_m)})$ are linearly independent over time for any initial condition in $\Omega$ or a parameter-dependent positively invariant subspace, and the mapping $r$ is injective in $\Theta$, the parameters are identifiable.

**Joint recovery of states and parameters.** If the above conditions hold, the pair $(x_0, \theta_0)$ can be reconstructed from the outputs, a finite number of their derivatives and a set of time instants which depends on the analytic nature of functions $\phi_{j,l} = g_{j,l}(\mathbf{y}^{(d'_1,...,d'_m)})$. In this case, we need to solve one of the following systems:

* Low-order derivatives, several time instants:

$$\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). \quad (13)$$

* High-order derivatives, one time instant:

$$\left( \begin{array}{ccc} \phi_{j,1}(\tilde{t}) & \cdots & \phi_{j,q_j}(\tilde{t}) \\ \frac{d\phi_{j,1}}{dt}(\tilde{t}) & \cdots & \frac{d\phi_{j,q_j}}{dt}(\tilde{t}) \\ \vdots & \ddots & \vdots \\ \frac{d^{q_j-1}\phi_{j,1}}{dt^{q_j-1}}(\tilde{t}) & \cdots & \frac{d^{q_j-1}\phi_{j,q_j}}{dt^{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}) \\ \frac{d\phi_{j,0}}{dt}(\tilde{t}) \\ \vdots \\ \frac{d^{q_j-1}\phi_{j,0}}{dt^{q_j-1}}(\tilde{t}) \end{array} \right). \quad (14)$$

Once we have the parameters, we can recover the state at one of the previously chosen times inverting the mapping $\mathcal{L}_{f_{\theta},h_{\theta},\{d_1,...,d_m\}}$.

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 joint 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 \ge 0$ (see, for example, [10, Chapter 10.4]). Then, we consider the following reduced system:

$$ \begin{cases} \dot{S} = -\beta SI + \mu(1 - S - I), \\ \dot{I} = \beta SI - \gamma I, \\ y = kI, \end{cases} \quad \forall t \ge 0, \quad (15) $$

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 \le 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 ([10, Chapter 10.4], [49]):

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

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

* The basic reproduction number is defined as $R_0 = \beta/\gamma$, which 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_f = (1, 0)^{\mathrm{T}}$, is globally asymptotically stable when $R_0 < 1$ and unstable when $R_0 > 1$; in particular, in this case, it is a saddle point whose stable manifold is $I \equiv 0$.

* When $R_0 > 1$, there exists an endemic equilibrium (EE) in $\Omega$ given by $P_e(\theta) = \left(\frac{1}{R_0}, \mu \frac{1 - 1/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 (15) are analytic in $\mathcal{I}$. We will see that the procedure described in Algorithm 2 is appropriate for this case, and gives more general results than Algorithm 1. Hence, in the following, we are going to focus on Algorithm 2 considering we know $y(x_0, \theta_0)$ in $[a, b) \subset \mathcal{I}$. In particular, to check the observability, identifiability or joint observability-identifiability of System (15), 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 (15) 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 (15) can be consider as a system of the form given in (2), with

$$f(x, \theta) = \begin{pmatrix} -\beta x_1 x_2 + \mu(1 - x_1 - x_2) \\ \beta x_1 x_2 - \gamma x_2 \end{pmatrix}, \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 (15):

$$\mathcal{L}_{\theta,d_1} : \xi \mapsto (y_{(\xi,\theta)}^{(0)}(0), \dots, y_{(\xi,\theta)}^{(d_1)}(0)).$$

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

$$y_{(\xi,\theta)}^{(0)}(t) = h_{\theta}(x(t; \xi, \theta)) = kI(t; \xi, \theta), \quad \forall t \ge 0.$$

Notice that $\mathcal{L}_{\theta,0}(\xi) = y_{(\xi,\theta)}^{(0)}(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 \to [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 (15) if, and only if, for any initial condition $\xi \in \Omega_1$, $S(t; \xi, \theta) = G(I(t; \xi, \theta))$, for all $t \ge 0$. Two examples of such orbits are $(S, I) \equiv P_f$ and, whenever $\mathcal{R}_0 > 1$, $(S, I) \equiv P_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 (15) that we will derive from it is less restrictive.

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

$$\begin{aligned} y_{(\xi,\theta)}^{(1)}(t) &= Dh_{\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_{(\xi,\theta)}^{(0)}(t), \quad \forall t \ge 0. \end{aligned}$$

Let us consider $\Omega_1 = \{(\xi_1, \xi_2)^T \in \Omega : \xi_2 \neq 0\} = \{(\xi_1, \xi_2)^T \in [0, 1] \times (0, 1] : \xi_1 + \xi_2 \le 1\}$.

**Lemma 5.** For any $\theta = (k, \beta, \gamma, \mu)^T \in \Theta$, the following map is injective:

$$\mathcal{L}_{\theta,1}: \Omega_1 \to \mathbb{R}^2 : \xi \mapsto (y_{(\xi,\theta)}^{(0)}(0), y_{(\xi,\theta)}^{(1)}(0)). \qquad (16)$$

*Proof.* The proof is straightforward, since, given $\theta = (k, \beta, \gamma, \mu) \in \Theta$, $\xi \in \Omega_1$,

$$y_{(\xi,\theta)}^{(0)}(0) = k\xi_2, \quad y_{(\xi,\theta)}^{(1)}(0) = k(\beta\xi_1 - \gamma)\xi_2.$$

Therefore, it is easy to see that (16) 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 2, we point out that $\Omega_1$ is positively invariant with respect to the system of ODEs of System (15). 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 (15) (recall that this is the stable manifold associated to $P_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 (15). Hence, System (15) 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 2.

**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 (15), and maps $g: \mathcal{D} \to \mathbb{R}^{q+p}$ and $r: \Theta \to \mathbb{R}^p$, for some suitable $\mathcal{D} \subset \mathbb{R}^{d_1+1}$, $d_1 \in \mathbb{N} \cup \{0\}$, $q, p \in \mathbb{N}$, such that (C1), (C2) and (C3) of Theorem 2 are satisfied for some connected $S \subset \mathcal{I}$ such that $[a, b) \subset S$ and $\{g_{j,l}(y_{(\xi,\theta)}^{(0)}(t), \dots, y_{(\xi,\theta)}^{(d_1)}(t))\}_{j=1,\dots,p}^{l=1,\dots,q_j}$ are analytic w.r.t. $t \in S$, for any $(\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 \ge 0, \qquad (17)$$

neither $d_1' = 0$ nor $d_1' = 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 (17).

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 \ge 0. \qquad (18)$$

Notice that, if $y \ne 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 (17) as

$$S = \frac{1}{\beta} \left( \frac{\dot{y}}{y} + \gamma \right), \quad I = \frac{y}{k}, \quad \forall t \ge 0. \qquad (19)$$

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

$$\ddot{y} - \frac{\dot{y}^2}{y} = -y\mu(\gamma - \beta) - y^2\frac{\beta}{k}(\gamma + \mu) - \dot{y}\mu - y\dot{y}\frac{\beta}{k}, \quad \forall t \ge 0. \qquad (20)$$

Then, we can see that (C1) of Theorem 2 is satisfied with $d_1' = 2$, $p = 1$, $q_1 = q = 4$,

$$g_0(y, \dot{y}, \ddot{y}) = \ddot{y} - \frac{\dot{y}^2}{y}, \quad g_1(y, \dot{y}, \ddot{y}) = -y, \quad g_2(y, \dot{y}, \ddot{y}) = -y^2, \quad g_3(y, \dot{y}, \ddot{y}) = -\dot{y}, \quad g_4(y, \dot{y}, \ddot{y}) = -y\dot{y},$$

and

$$r(\theta) = \left( \mu(\gamma - \beta), \frac{\beta}{k}(\gamma + \mu), \mu, \frac{\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 \to r(\Theta) \subset \mathbb{R}^4$ is injective. Indeed, given some parameters vector $\theta = (k, \beta, \gamma, \mu)^T \in \Theta$, it is easy to see that we can invert the equation $r(\theta) = \sigma$ and obtain uniquely

$$(k, \beta, \gamma, \mu) = \left( \frac{1}{\sigma_4} \left( \frac{\sigma_2}{\sigma_4} - \sigma_3 - \frac{\sigma_1}{\sigma_3} \right), \frac{\sigma_2}{\sigma_4} - \sigma_3 - \frac{\sigma_1}{\sigma_3}, \frac{\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 $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 (15). 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_e(\theta)\}) \times \{\theta\}$ and $(\xi, \theta) \in \Gamma_{2,\Theta}$, i.e., $\Omega_{2,\theta} = \Omega_1 \setminus \{P_e(\theta)\}$, where $P_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_1y(t) + a_2y(t)^2 + a_3\dot{y}(t) + a_4y(t)\dot{y}(t) = 0, \quad \forall t \in [a, b). \quad (21)$$

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 (21) is true for values $\{a_i\}_{i=1}^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 (15) with initial condition $\xi$ and parameters vector $\theta$, recall that $\Omega_1$ is positively invariant w.r.t. the ODE system given in (15), i.e., $(S(t), I(t))^{\mathrm{T}} \in \Omega_1$, for all $t \ge 0$. Then, we can rewrite (21) as

$$\begin{gathered} a_1kI(t) + a_2k^2I^2(t) + a_3kI(t)(\beta S(t) - \gamma) + a_4k^2I^2(t)(\beta S(t) - \gamma) = 0 \implies \\ (S, I)^{\mathrm{T}} \in \Omega_1 \implies a_1 + a_2kI(t) + a_3(\beta S(t) - \gamma) + a_4kI(t)(\beta S(t) - \gamma) = 0, \quad \forall t \in [a, b), \end{gathered}$$

which leads to

$$(a_1 - a_3\gamma) + a_3\beta S(t) + (a_2k - a_4k\gamma)I(t) + a_4k\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), \quad (22)$$

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, (22) 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_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_e(\theta)$; on the other hand, if $\dot{S} \equiv 0$, then, $S \equiv \xi_1$ and

$$I \equiv \frac{\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_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 (22) 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) = -\frac{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) \frac{BC - AD}{(C + DS(t))^2} = \left( (\beta S(t) + \mu) \frac{A + BS(t)}{C + DS(t)} - \mu(1 - S(t)) \right) \frac{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)) \frac{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) \frac{A + BS(t)}{C + DS(t)} - \mu(1 - S(t)) \right) \frac{BC - AD}{(C + DS(t))^2} = (\gamma - \beta S(t)) \frac{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

$$ c_0 = \mu ABC - \mu A^2 D + \mu ACD - \gamma AC^2 - \mu BC^2, $$

$$ c_1 = \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, $$

$$ c_2 = \beta BC^2 - (2\gamma - \mu) BCD - (\gamma + \mu) AD^2 - \beta ABD + 2\beta ACD + \beta B^2 C, $$

$$ c_3 = \beta AD^2 - \gamma BD^2 + 2\beta BCD, $$

$$ c_4 = \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 (22) 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} = \{ (\xi_1, \xi_2)^T \in [0, 1) \times (0, 1] : \xi_1 + \xi_2 \le 1 \} \setminus \{ P_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 (15).
**Lemma 7.** Let $\theta \in \Theta$. The set $\Omega_{2,\theta} = \{(\xi_1, \xi_2)^T \in [0,1) \times (0,1] : \xi_1 + \xi_2 \le 1\} \setminus \{P_e(\theta)\}$ is positively invariant with respect to the system

$$ \left\{ \begin{array}{l} \dot{S} = -\beta SI + \mu(1 - S - I), \\ \dot{I} = \beta SI - \gamma I. \end{array} \right. \qquad (23) $$

*Proof.* As said at the beginning of Section 4 the set

$$ \Omega = \{(\xi_1, \xi_2)^T \in [0,1]^2 : \xi_1 + \xi_2 \le 1\} $$

is positively invariant with respect to System (23).

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

On the other hand, if an initial condition $\xi$ is such that $\xi_2 \ne 0$, then $I(t) > 0$, for all $t \ge 0$, due to the uniqueness of solutions. Indeed, since $S(t) \ge 0$, $\forall t \ge 0$,

$$ \dot{I}(t) = (\beta S(t) - \gamma)I(t) \ge -\gamma I(t), \quad \forall t \ge 0 \quad \implies \quad I(t) \ge \xi_2 e^{-\gamma t} > 0, \quad \forall t \ge 0. $$

Thus, $\{(\xi_1, \xi_2)^T \in [0,1) \times (0,1] : \xi_1 + \xi_2 \le 1\} \setminus \{P_e(\theta)\}$ is also a positively invariant set with respect to System (23). □

Hence, Step 2 is finished, which allows us to conclude that System (15) 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 2.

**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_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 \text{constant} \ne 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 (15) 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 2 or Steps 4–9 in Algorithm 3.

**Remark 8.** Applying the DAISY algorithm [5] to this model yields the following differential relations between the states, the output, and the parameters:

$$ 0 = y - kI, $$

$$ 0 = \dot{y} - \beta Sy + \gamma y, $$

$$ 0 = kyy - k\dot{y}^2 + \beta y^2 \dot{y} + k\mu y\dot{y} + \beta(\gamma + \mu)y^3 + k\mu(\gamma - \beta)y^2. $$

These equations are equivalent to (19) and (20) when $y \ne 0$. The importance now is finding a suitable rewriting of the last equation such that (C2) and (C3) of Theorem 2 can be checked. In particular, it is clear that all functions in this equation are linearly dependent, so (C3) is not satisfied if we consider the equation written this way. This moreover highlights that (20) may not be the unique suitable rewriting. This example illustrates how our methodology complements symbolic approaches such as DAISY.
## 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 (15), 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:

$$ egin{cases} \dot{S} = -\beta SI, \\ \dot{I} = \beta SI - \gamma I, \\ y = kI, \end{cases} \qquad (24) $$

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 (24) is observable, identifiable or jointly observable-identifiable, it seems it will satisfy the same properties as the previous System (15), 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 (24) is observable on $\Omega_1 = \{(\xi_1, \xi_2)^T \in [0, 1] \times (0, 1] : \xi_1 + \xi_2 \le 1\}$ in any $[a, b) \subset I$ with parameters in $\Theta = (0, 1] \times (0, \infty)^2$. Regarding Step 2, as in (19) and (20) with $\mu = 0$, we obtain that

$$ I = \frac{y}{k}, \quad S = \frac{1}{\beta} \left( \frac{\dot{y}}{y} + \gamma \right) \qquad (25) $$

and

$$ \ddot{y} - \frac{\dot{y}}{y} = -\frac{\beta\gamma}{k}y^2 - \frac{\beta}{k}y\dot{y}. \qquad (26) $$

Then,

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

which is clearly not injective in $\Theta$, and hence we cannot complete Step 2 using (26). However, recalling Remark 1, System (24) 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 (24) is identifiable on $\Theta = (0, 1] \times (0, \infty)^2$ in any $[a, b) \subset I$, with initial conditions in $\tilde{\Omega} = \{(\xi_1, \xi_2)^T \in (0, 1)^2 : \xi_1 + \xi_2 \le 1\}^2$.

*Proof.* System (24) 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$ uniquely.

First of all, one can prove, similarly to Lemma 7, that $\tilde{\Omega}$ is positively invariant with respect to the ODE system in System (24). Equations (25) and (26) 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} - \frac{\dot{y}}{y} = -\sigma_1 y^2 - \sigma_2 y \dot{y}, \qquad (27) $$

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 \neq 0$ in $[a, b)$, and hence $y \neq 0$, which implies

$$ a_1 y + a_2 \dot{y} = 0. $$

<footnote n="2">²Note that $\tilde{\Omega}$ is independent of the choice of $\theta$.</footnote>
This is equivalent to

$$a_1 kI + a_2 kI(\beta S - \gamma) = 0, \forall t \in [a, b), \text{ which implies } 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 \ge 0$, which cannot occur in $\tilde{\Omega}$. Hence, $y^2$ and $yy$ are linearly independent when we consider $\xi \in \tilde{\Omega}$ and $\theta \in \Theta$, and therefore there is a unique solution $\sigma$ to (27). Then, we can determine uniquely $\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 = \frac{y(0)}{\xi_2}, \quad \gamma = \frac{\sigma_1}{\sigma_2}, \quad \beta = \frac{1}{\xi_1} \left( \frac{\dot{y}(0)}{y(0)} + \frac{\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 (24): $X = \beta S$ and $Y = kI$. Then,

$$\begin{cases} \dot{X} = -\sigma_2 XY, \\ \dot{Y} = Y(X - \gamma), \end{cases}$$

and we can recover $\beta\xi_1$ and $k\xi_2$ integrating backwards considering as initial condition

$$X(a) = \frac{\dot{y}(a)}{y(a)} + \frac{\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 uniquely $k, \beta$ and $\gamma$ as

$$k = \frac{Y(0)}{\xi_2}, \quad \beta = \frac{X(0)}{\xi_1}, \quad \gamma = \frac{\sigma_1}{\sigma_2}.$$

□

Therefore, System (24) 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 I$. Nevertheless, it is not jointly observable-identifiable on $\tilde{\Omega} \times \Theta$. Actually, in [19], 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. This can be as well checked with DAISY and StructuralIdentifiability.jl softwares. In fact, considering (25) and (26), if $(\sigma_1, \sigma_2)$ is the solution to

$$\ddot{y} - \frac{\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 = \frac{y}{k_0}, \quad S = \frac{1}{\beta_0} \left( \frac{\dot{y}}{y} + \gamma_0 \right), \quad \frac{\beta_0 \gamma_0}{k_0} = \sigma_1, \quad \frac{\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$ uniquely as follows:

$$\gamma_0 = \frac{\sigma_1}{\sigma_2}, \quad \frac{\beta_0}{k_0} = \sigma_2, \quad \beta_0 S_0 = \frac{\dot{y}(0)}{y(0)} + \frac{\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)^T$ and $(k_2, \beta_2, \gamma_2, \mu_2)^T$, $\mu_2 \neq 0$. To tackle this question, we consider an *extended SIRS model* with parameters

$$ \theta := (k, \beta, \gamma, \mu)^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 \to 0$.

Then again, in Step 2 we would obtain

$$ r(k, \beta, \gamma, \mu) = \left( \mu(\gamma - \beta), \frac{\beta}{k}(\gamma + \mu), \mu, \frac{\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_{\text{SIR}}(t)$ and $x_{\text{SIRS}}(t)$ be the solutions to the SIR given by the ODE system of (24) and the SIRS model (23), respectively, with the same initial condition, and we study the dependence of $||x_{\text{SIR}}(t) - x_{\text{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 [49], particularized to our autonomous context:

**Theorem 4.** [49, 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), \quad y(t_0) = y_0, \quad \text{and} \quad \dot{z} = f(z) + g(z), \quad 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)) + \frac{\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_{\text{SIR}}(t) - x_{\text{SIRS}}(t)||_2$.

Notice that both models share the same positively invariant set

$$ \Omega = \{(\xi_1, \xi_2)^T \in [0, 1]^2 : \xi_1 + \xi_2 \leq 1\}, $$

which is compact, but the derivatives are well defined. Let $x_{\text{SIR}} = (S_{\text{SIR}}, I_{\text{SIR}})^T$ and $x_{\text{SIRS}} = (S_{\text{SIRS}}, I_{\text{SIRS}})^T$ be the state variables associated to the ODE system of (24) and (23), 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)^T$ and $\nu = (\beta, \gamma)^T$ two parameters vectors, and the following function:

$$ f_{\text{SIR}}(x, \nu) = \begin{pmatrix} -\beta x_1 x_2 \\ \beta x_1 x_2 - \gamma x_2 \end{pmatrix}, \quad x \in \Omega. $$
Notice that $f_{\text{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) = \begin{pmatrix} \mu(1 - x_1 - x_2) \\ 0 \end{pmatrix}, \quad x \in \Omega.$$

Then, $g$ satisfies

$$\|g(x, \mu)\|_2 = \left\| \begin{pmatrix} \mu(1 - x_1 - x_2) \\ 0 \end{pmatrix} \right\|_2 \le \mu, \quad \forall x \in \Omega.$$

We consider now the two following systems:

$$\begin{pmatrix} \dot{S}_{\text{SIR}}(t; x_0, \nu) \\ \dot{I}_{\text{SIR}}(t; x_0, \nu) \end{pmatrix} = f_{\text{SIR}}(x_{\text{SIR}}(t; x_0, \nu), \nu)$$

and

$$\begin{pmatrix} \dot{S}_{\text{SIRS}}(t; x_0, \eta) \\ \dot{I}_{\text{SIRS}}(t; x_0, \eta) \end{pmatrix} = f_{\text{SIR}}(x_{\text{SIRS}}(t; x_0, \eta), \nu) + g(x_{\text{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_{\text{SIR}}(t) = x_{\text{SIR}}(t; x_0, \nu)$ and $x_{\text{SIRS}}(t) = x_{\text{SIR}}(t; x_0, \eta)$. It is fulfilled that

$$\|x_{\text{SIR}}(t) - x_{\text{SIRS}}(t)\|_2 \le \frac{\mu}{L}[\exp(L(t - t_0)) - 1], \quad \forall t \ge t_0.$$

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

$$\|x_{\text{SIR}}(t) - x_{\text{SIRS}}(t)\|_2 \le \frac{\mu}{L}[\exp(L(t - t_0)) - 1] < \varepsilon, \quad \forall t \in [t_0, \bar{t}]. \qquad (28)$$

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_{\text{SIR}}(0) = x_{\text{SIRS}}(0) = (0.9, 0.1)^{\text{T}}$. Moreover, the infectious compartment $I_{\text{SIRS}}$ presents a slow-fast behavior for the SIRS model when near to the invariant manifold $I_{\text{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_{\text{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.

### 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}, \dot{y}^2\}$ 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)^{\text{T}} \in (0, 1)^2 : \xi_1 + \xi_2 \le 1\}$. Then, given some observations $y(x_0, \theta_0)$, we can perform Step 3 of Algorithm 2: 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 (28), such that there exists a unique solution $\sigma$
<figure><img src="image_2.png" /><figcaption>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).</figcaption></figure>

fulfilling

$$ \begin{pmatrix} 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{pmatrix} \begin{pmatrix} \sigma_1 \\ \sigma_2 \\ \sigma_3 \\ \sigma_4 \end{pmatrix} = \begin{pmatrix} \frac{\dot{y}^2(t_1)}{y(t_1)} - \ddot{y}(t_1) \\ \frac{\dot{y}^2(t_2)}{y(t_2)} - \ddot{y}(t_2) \\ \frac{\dot{y}^2(t_3)}{y(t_3)} - \ddot{y}(t_3) \\ \frac{\dot{y}^2(t_4)}{y(t_4)} - \ddot{y}(t_4) \end{pmatrix} $$

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 (26), i.e.,

$$ \ddot{y} - \frac{\dot{y}}{y} = -\frac{\beta\gamma}{k}y^2 - \frac{\beta}{k}y\dot{y}, \quad \text{which may be rewritten as} \quad \frac{d}{dt}\frac{\dot{y}}{y} + \frac{\beta\gamma}{k}y + \frac{\beta}{k}\dot{y} = 0, $$

i.e., $\left\{ \frac{d}{dt}\frac{\dot{y}}{y}, y, \dot{y} \right\}$ are linearly dependent for the SIR model whenever $I_{SIR} \neq 0$. However, this is not true for the SIRS model. Indeed, consider System (15). Let $a_1, a_2, a_3 \in \mathbb{R}$ such that

$$ a_1 \frac{d}{dt} \frac{\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_2kI + a_3kI(\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}, \frac{d\dot{y}}{dt} \frac{d\dot{y}}{dy}$ may help us determine the model.

These two approaches highlight the importance of the structural equations associated to our models, which can be exploited for model discrimination in particular cases.

# 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 case, we show how joint observability-identifiability can be established using our proposed framework under different modeling assumptions and data limitations. For clarity, we systematically omit the equation for the Recovered compartment when the total population is constant, as it does not affect the theoretical results. To maintain focus on the theoretical aspects, detailed computational steps are omitted but can be derived directly following the procedures outlined in previous sections.

It is worth noting that the following examples can also be analyzed using symbolic tools developed for rational systems, such as DAISY [5] and StructuralIdentifiability.jl softwares [25], such as performed, for instance, in [17, 55]. However, our interest here lies in explicitly determining the sets of initial conditions for which joint observability-identifiability holds—information that these algorithms do not provide.

## The SIR model with demography

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

$$ egin{cases} \dot{S} = \delta N - \beta \frac{SI}{N} - \delta S, \\ \dot{I} = \beta \frac{SI}{N} - (\gamma + \delta)I, & \forall t \ge 0, \\ y = kI, \end{cases} $$

where $\delta > 0$ denotes the birth and death rates, assumed equal, $N$ is the total population (constant), and $\beta, \gamma, 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:

$$ egin{cases} \dot{S} = \delta - \beta SI - \delta S, \\ \dot{I} = \beta SI - (\gamma + \delta)I, & \forall t \ge 0, \\ y = kI, \end{cases} \quad (29) $$

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

For this normalized system, we set $\Omega = \{(\xi_1, \xi_2)^T \in [0, 1]^2 : \xi_1 + \xi_2 \le 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 \Rightarrow \quad S = \frac{1}{\beta} \left( \frac{\dot{y}}{y} + \gamma + \delta \right), \quad I = \frac{y}{k}, $$
assuming $y \neq 0$, and

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

Then, we consider $d_1 = 1$, $d'_1 = 2$, $p = 1$, $q_1 = q = 4$,

$$ g_0(y, \dot{y}, \ddot{y}) = \ddot{y} - \frac{\dot{y}^2}{y}, \quad g_1(y, \dot{y}, \ddot{y}) = y, \quad g_2(y, \dot{y}, \ddot{y}) = -y^2, \quad g_3(y, \dot{y}, \ddot{y}) = -\dot{y}, \quad g_4(y, \dot{y}, \ddot{y}) = -y\dot{y} $$

and

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

With these expressions, we confirm that the assumptions in Step 1 and Step 2 of Algorithm 2 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 \le 1\} \setminus \{P_e(\theta)\}$, for each $\theta \in \Theta$, where $P_e(\theta)$ denotes the equilibrium point associated with $\theta$, and $\Omega_{\theta}$ is positively invariant with respect to the ODE system of (29) when we consider $\theta$ as parameters vector.

Therefore, we conclude that System (29) 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}{l} \dot{S} = -\beta SI - \nu S, \\ \dot{I} = \beta SI - \gamma I, \\ \dot{V} = \nu S, \\ y = \nu(1 - V), \end{array} \quad \forall t \ge 0, \right. \qquad (30) $$

for some initial condition $(S(0), I(0), V(0))^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)^T \in [0, 1]^3 : \xi_1 + \xi_2 + \xi_3 \le 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 = -\frac{\dot{y}}{\nu^2}, \quad I = -\frac{1}{\beta}\left(\frac{\ddot{y}}{\dot{y}} + \nu\right), \quad V = 1 - \frac{y}{\nu}, $$

assuming $\dot{y} \neq 0$, and

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

Then, we consider $d_1 = 2$, $d'_1 = 3$, $p = 1$, $q_1 = q = 4$,

$$ g_0(y, \dot{y}, \ddot{y}, y^{(3)}) = y^{(3)} - \frac{\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, \frac{\beta}{\nu}, \gamma, \frac{\beta}{\nu^2} \right).$$

Following the procedures established in earlier sections, we confirm that the assumptions of Step 1 and Step 2 in Algorithm 2 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 \le 1\}$, where $\Omega$ does not depend on $\theta \in \Theta$ and is positively invariant under the ODE system defined by (30).

Thus, we conclude that System (30) 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 (24)) 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{aligned} \dot{S} &= -\beta SI, \\ \dot{I} &= \beta SI - \gamma I, \\ y &= \beta SI, \end{aligned} \right. \qquad (31)$$

for some initial condition $(S(0), I(0)) \in \Omega = \{(\xi_1, \xi_2)^T \in [0, 1]^2 : \xi_1 + \xi_2 \le 1\}$, which is positively invariant with respect to the system of ODEs given in (31). 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{aligned} y &= \beta SI, \\ \dot{y} &= -y\beta I + y\beta S - \gamma y. \end{aligned} \right.$$

Notice that this system is nonlinear in $S$ and $I$, making it challenging to directly apply Step 1 of Algorithm 2. 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{aligned} \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{aligned} \right. \quad \Longrightarrow \quad \left\{ \begin{aligned} S &= \frac{1}{\beta} \frac{\dot{y}}{y} + \frac{1}{\beta\gamma} \left( \frac{\ddot{y}}{y} - \frac{\dot{y}^2}{y^2} \right) + \frac{2}{\gamma} y + \frac{\gamma}{\beta}, \\ I &= \frac{1}{\beta\gamma} \left( \frac{\ddot{y}}{y} - \frac{\dot{y}^2}{y^2} \right) + \frac{2}{\gamma} y. \end{aligned} \right.$$

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

$$- \left( \frac{d}{dt} \frac{\dot{y}}{y} \right)^2 = \gamma \frac{\dot{y}}{y} \frac{d}{dt} \frac{\dot{y}}{y} + 2\beta\gamma\dot{y} + \beta\gamma^2 y + 4\beta y \frac{d}{dt} \frac{\dot{y}}{y} + 4\beta^2 y^2 + \gamma^2 \frac{d}{dt} \frac{\dot{y}}{y}.$$

Then, we set $d_1 = d'_1 = 2$, $p = 1$, $q_1 = q = 6$,

$$g_0(y, \dot{y}, \ddot{y}) = - \left( \frac{d}{dt} \frac{\dot{y}}{y} \right)^2, \quad g_1(y, \dot{y}, \ddot{y}) = \frac{\dot{y}}{y} \frac{d}{dt} \frac{\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 \frac{d \dot{y}}{dt y}, \quad g_5(y, \dot{y}, \ddot{y}) = y^2, \quad g_6(y, \dot{y}, \ddot{y}) = \frac{d \dot{y}}{dt y},$$

and

$$r(\theta) = (\gamma, 2\beta\gamma, \beta\gamma^2, 4\beta, 4\beta^2, \gamma^2).$$

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

Consequently, we conclude that System (31) 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-Infected-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 (30). 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}{l}
\dot{S} = A - \beta SI - (\nu + \delta)S, \\
\dot{I} = \beta SI - \delta I, \\
\dot{V} = \nu S - \delta V, \\
y = \begin{pmatrix} y_1 \\ y_2 \end{pmatrix} = \begin{pmatrix} \nu(S + I) \\ \delta(S + I + V) \end{pmatrix},
\end{array}
\right.
\quad \forall t \ge 0, \qquad (32)$$

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 (30).

The natural positively invariant set for this system is $\Omega = \{(\xi_1, \xi_2, \xi_3)^T \in [0, A/\delta]^3 : \xi_1 + \xi_2 + \xi_3 \le 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 (30). 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}{l} 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}{l} S = \frac{A}{\nu} - \frac{\delta}{\nu^2} y_1 - \frac{\dot{y}_1}{\nu^2}, \\ I = \left( \frac{\nu + \delta}{\nu^2} \right) y_1 - \frac{A}{\nu} + \frac{\dot{y}_1}{\nu^2}, \\ V = \frac{y_2}{\delta} - \frac{y_1}{\nu}, \end{array} \right. $$

and

$$ \left\{ \begin{array}{l} \dot{y}_1^2 = A\nu^2 \frac{\delta\nu - A\beta}{\beta} + \nu \left( A\nu + 2A\delta - \frac{\delta\nu^2 + \delta^2\nu}{\beta} \right) y_1 + \nu \left( 2A - \frac{\nu^2 + 2\delta\nu}{\beta} \right) \dot{y}_1 \\ \qquad -\delta(\nu + \delta)y_1^2 - (\nu + 2\delta)y_1\dot{y}_1 - \frac{\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{aligned} \dot{y}_1^2 + \sigma_2^2 y_1^2 - 2y_1 \dot{y}_1 \sigma_2 &= A\nu^2 \frac{\delta\nu - A\beta}{\beta} + \nu \left( A\nu + 2A\delta - \frac{\delta\nu^2 + \delta^2\nu}{\beta} \right) y_1 \\ &\quad + \nu \left( 2A - \frac{\nu^2 + 2\delta\nu}{\beta} \right) \dot{y}_1 - \nu y_1 (-\sigma_2 y_1 + \dot{y}_1) - \frac{\nu^2}{\beta} \ddot{y}_1, \\ \dot{y}_2 &= A\delta - \delta y_2. \end{aligned} \right. $$

This reformulation ensures the linear independence condition in Step 2 of Algorithm 2. Then, we consider $d_1 = 1, d_2 = 0, d'_1 = 2, d'_2 = 1, p = 2, q_1 = 5, q_2 = 2,$

$$ g_{1,0}(y_1, \dot{y}_1, \ddot{y}_1, y_2, \dot{y}_2) = \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, $$

$$ g_{1,2}(y_1, \dot{y}_1, \ddot{y}_1, y_2, \dot{y}_2) = y_1, \quad g_{1,3}(y_1, \dot{y}_1, \ddot{y}_1, y_2, \dot{y}_2) = \dot{y}_1, $$

$$ g_{1,4}(y_1, \dot{y}_1, \ddot{y}_1, y_2, \dot{y}_2) = -y_1(-\sigma_2 y_1 + \dot{y}_1), \quad g_{1,5}(y_1, \dot{y}_1, \ddot{y}_1, y_2, \dot{y}_2) = -\ddot{y}_1, $$

$$ g_{2,0}(y_1, \dot{y}_1, \ddot{y}_1, y_2, \dot{y}_2) = \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 \frac{\delta\nu - A\beta}{\beta}, \nu \left( A\nu + 2A\delta - \frac{\delta\nu^2 + \delta^2\nu}{\beta} \right), \nu \left( 2A - \frac{\nu^2 + 2\delta\nu}{\beta} \right), \nu \frac{\nu^2}{\beta}, A\delta, \delta \right). $$

Given this, the assumptions in Step 1 and Step 2 of Algorithm 2 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_e(\theta)\}$, where $P_e(\theta)$ is the equilibrium point of the system and $\Omega_\theta$ is positively invariant under the ODEs of System (32) when considering $\theta$ as the parameters vector.

Therefore, we conclude that System (32) 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 and practical behavior of the proposed framework given in Sections 4 and 3 for observability, identifiability, and joint observability.
ability-identifiability in epidemiological models. The simulations serve both to confirm the theoretical results established in Sections 3 and 4 and to provide intuition on the numerical conditioning of the recovery procedure. In particular, they highlight how the spacing of the selected observation times affects the stability of the linear systems used for parameter reconstruction, and how the estimated parameters converge to their theoretical values when the model is identifiable. Two representative examples are analyzed, one based on an SIRS model and another on an SIR model, each chosen to exhibit distinct epidemiological dynamics and to validate the results in Section 5:

* **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 $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, $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_0I$, 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 3 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., [80] and [81]). 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 3, omitting that of Algorithm 2 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 the *continuous* observation of $y = k_0I$ 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_{\text{vec}} = (t_i)_{i \in \{0, \dots, N\}}$, where $t_i = ih$ and $N = T_{\text{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} = \frac{\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 3 to estimate the system parameters and initial condition. This method requires selecting a time $\tilde{t} \in t_{\text{vec}}$ such that the following system has a unique solution (see (10)-(11)):

$$\left( \frac{d^k y}{dt^k} \quad \frac{d^k y^2}{dt^k} \quad \frac{d^k \dot{y}}{dt^k} \quad \frac{d^k y \dot{y}}{dt^k} \right)_{\substack{t=\tilde{t} \\ k=0,1,2,3}} \sigma^T = \left( \frac{d^k}{dt^k} \left( \frac{\dot{y}^2}{y} - \ddot{y} \right) \right)_{\substack{t=\tilde{t} \\ k=0,1,2,3}} . \quad (33)$$

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), \frac{\beta_0}{k_0}(\gamma_0 + \mu_0), \mu_0, \frac{\beta_0}{k_0} \right).$$

The initial conditions $S_0$ and $I_0$ are then computed as

$$S_0 = \frac{1}{\beta_0} \left( \frac{\dot{y}(0)}{y(0)} + \gamma_0 \right) \quad \text{and} \quad I_0 = \frac{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.83) \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_{\text{max}}/h+1 = 161$. For each time point in $t_{\text{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 (33), which should be non-zero. The values obtained were of the order of $10^{-17}$.

* The condition number of the matrix in (33) (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 parameters 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_{OLS} = (k_{OLS}, \beta_{OLS}, \gamma_{OLS}, \mu_{OLS}, S_{0,OLS})$ by solving

$$ \theta_{OLS} = \underset{\theta \in \Theta}{\operatorname{argmin}} \sum_{i=0}^{T_{\max}} A (\theta_1 I(t_i; \theta) - y(t_i))^2, \quad t_i = i, $$

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

$$ \begin{cases} \dot{S}(t; \theta) = -\theta_2 S(t; \theta) I(t; \theta) + \theta_4 (1 - S(t; \theta) - I(t; \theta)), \\ \dot{I}(t; \theta) = \theta_2 S(t; \theta) I(t; \theta) - \theta_3 I(t; \theta), \end{cases} $$

with $(S(0; \theta), I(0; \theta)) = \left(\theta_5, \frac{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 [60]). 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 $\mathbb{R}_0$ estimates for diseases such as smallpox, pertussis, or COVID-19 (e.g., [31], [51], [54], [53], [78], [79]).

We initialize the OLS algorithm several times with random initial conditions, denoted as I.C., drawn from a uniform distribution $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 U(0, 1)$. For conciseness, we report results from 5 representative tests. Table 1 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 parameters vector $\theta_0$), Obj. value (final objective function value), and Time (s) (computational time for each test in seconds).
<table><caption>Table 1: Results for Scenario 2 applied to Case 2 using the OLS method.</caption><thead><tr><th>Test</th><th>I.C.</th><th>Abs. error θ<sub>0</sub></th><th>Obj. value</th><th>Time (s)</th></tr></thead><tbody><tr><td>1</td><td>(0.545 1.967 0.414 0.82 0.718)</td><td>(0.019 0.015 6.8e-6 4.7e-6 0.052)</td><td>6.324e-12</td><td>116.517</td></tr><tr><td>2</td><td>(0.97 1.599 0.332 0.106 0.611)</td><td>(0.226 0.189 8.9e-12 9.3e-13 0.387)</td><td>7.657e-19</td><td>1.904</td></tr><tr><td>3</td><td>(0.785 1.276 0.1 0.266 0.154)</td><td>(0.25 0.208 4e-11 3.6e-12 0.409)</td><td>9.036e-18</td><td>1.546</td></tr><tr><td>4</td><td>(0.686 0.874 0.675 0.695 0.068)</td><td>(0.7 0.286 0.270 0.096 0.024)</td><td>0.003</td><td>621.808</td></tr><tr><td>5</td><td>(0.277 0.68 0.671 0.844 0.344)</td><td>(0.7 0.287 0.271 0.097 0.024)</td><td>0.006</td><td>1380.792</td></tr></tbody></table>

<table><caption>Table 2: Results for Scenario 2 applied to Case 2 using the OLS method: Approximations of the parameters.</caption><thead><tr><th>Test</th><th>Approx. of θ<sub>0</sub></th><th>Approx. of θ̃<sub>0</sub></th></tr></thead><tbody><tr><td>1</td><td>(0.319 0.265 0.1 4.7e-6 0.848)</td><td>(0.1 0.833 0.225)</td></tr><tr><td>2</td><td>(0.526 0.439 0.1 9.3e-13 0.513)</td><td>(0.1 0.833 0.225)</td></tr><tr><td>3</td><td>(0.55 0.458 0.1 3.6e-12 0.491)</td><td>(0.1 0.833 0.225)</td></tr><tr><td>4</td><td>(1 0.536 0.370 0.096 0.924)</td><td>(0.370 0.536 0.495)</td></tr><tr><td>5</td><td>(1 0.537 0.371 0.097 0.924)</td><td>(0.371 0.537 0.496)</td></tr></tbody></table>

The target parameters 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.83, 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 joint identifiability-observability in a broad class of systems of ODEs encompassing many classical epidemiological models, under the assumption of continuous, noise-free observations. Several symbolic differential-algebraic approaches, such as DAISY [5]
and StructuralIdentifiability.jl [25], have provided rigorous algebraic tools to determine structural identifiability in rational dynamical systems. For nonrational, analytic systems, we can find the software STRIKE-GOLDD [75], based on a generalization of the Observability Rank Condition. However, all these methods yield qualitative verdicts only, without identifying the regions of the parameters space or the sets of initial conditions where full recoverability holds. Building upon these advances, our study introduces a complementary and constructive framework that extends the analysis to general nonlinear (not necessarily analytic) systems and explicitly characterizes the conditions under which both parameters and initial states can be uniquely recovered.

This problem was previously proposed in [52], where the foundation of the approach was developed in a general context. In Section 3, we extended that framework by considering parameter-dependent sets of initial conditions and presenting Algorithms 1–3, which provide constructive procedures to recover the unknown parameters and/or initial conditions. These methods rely on solving a series of linear systems at suitably chosen time instants, offering a practical and systematic approach to recovering model parameters. If observations are available at time 0, the initial condition can be determined directly; otherwise, it can be reconstructed by integrating backwards using the estimated parameters. Furthermore, our approach can be naturally extended to models with piecewise constant parameters by allowing parameter updates at predetermined time instants.

One of the core ideas underlying this theoretical framework is the analysis of linear independence among some sets of functions, a concept seldom emphasized in the literature, where existing results typically address more restrictive settings (e.g., rational systems or models with known initial data). This approach has proven to be essential for establishing joint observability-identifiability, while reducing the dependence on high-order derivatives. As a consequence, it potentially enables more robust algorithms in the presence of noise compared to classical conditions such as the Hermann–Krener rank criterion [39].

In Section 4, we applied this framework to the SIRS model (with $\mu > 0$) when only a fraction of infected individuals is observed and proved that the system is jointly observable-identifiable under idealized observations. Specifically, we showed that observing only a fraction of the infected population suffices to determine all model parameters ($k, \beta, \gamma$, and $\mu$) and the initial condition of both susceptible and infected individuals. It is demonstrated that considering parameter-dependent sets of initial conditions is crucial in some cases. In contrast, for the limiting SIR model ($\mu = 0$), we demonstrated that the system remains observable and identifiable but loses joint observability-identifiability, since only certain parameter combinations ($\gamma, \beta/k, \beta S_0, kI_0$) can be recovered. This highlights the critical role of the SIRS model's cyclic dynamics (loss of immunity) in enabling full recoverability. Based on these findings, we proposed a theoretical basis for discriminating between the SIR and SIRS structures using short-time data.

In Section 5, we presented additional examples to illustrate the versatility of our approach across different model architectures and observation schemes. In these examples, we explicitly characterize the sets of initial conditions where joint observability-identifiability holds, a feature not provided by symbolic approaches, complementing these existing methodologies.

In Section 6, we validated the theoretical results through numerical experiments on both SIR and SIRS models under two different observational scenarios:

* **Scenario 1: Continuous** observations, assuming perfect knowledge of derivatives, where Algorithm 3 successfully approximated the parameter values in all SIRS cases.

* **Scenario 2:** Discrete daily observations, where Ordinary Least Squares was used for calibration in the SIR case. The identifiable parameter combinations ($\gamma, \beta/k, \beta S(0)$) and $\mu \approx 0$ were consistently captured in most of the tests. Occasional fits resembling the SIRS structure further highlighted the practical challenges of discriminating between SIR and SIRS dynamics during early epidemic stages.
This scenario illustrates the results in Section 4.2.

**Perspectives and possible extensions.** The proposed framework provides a constructive foundation that can be explored to be further extended to more complex epidemic models featuring time-varying parameters, delays, stochastic effects, or spatial heterogeneity. In particular, an effective strategy could be transforming these more complex systems into a system that falls within our framework.

We commented in Remark 7 the possibility of extending this framework to piecewise constant parameters under the knowledge of the time instants when their values change. Further, if parameters evolve continuously, they can be regarded as additional dynamic states. In particular, if their dynamics are governed by slow or piecewise constant processes, the model may preserve the structure of the presented framework.

On the other hand, systems of partial differential equations are also widely used for epidemiological modeling to consider, e.g., spatial heterogeneity or age of individuals. These models can be transformed into or approximated by systems of ODEs using different techniques ([18, 21, 30]).

Third, considering stochastic effects is of paramount importance due to the imperfect nature of data in physical and biological processes. In [12], authors propose to study the structural identifiability of the noise-free model along with an ODE system of the statistical moments of the involved stochastic processes, which again recovers the ODE-based framework.

Finally, for systems featuring delays, either constant or distributed (modeled through integro-differential equations), one can make use of the so-called *Linear Chain Trick* to approximate these delayed systems to a system of ODEs. This is performed by the inclusion of linear chains of new states of the system, leading to Erlang distributed stay times (see [14, 57, 65]).

Future work should as well explore extensions to scenarios with more realistic features, such as bounded measurement ranges (e.g., related to *interval observers*, see [68]).

## Acknowledgments

Authors ABK, BI and AMR were supported by the Spanish Government under projects PID2019-106337GB-I00, PID2023-146754NB-I00, and the European M-ERA.Net under project PCI2024-153478. Author AR was supported by the French National Research Agency (ANR) under project ANR NOCIME (New Observation and Control Issues Motivated by Epidemiology) 2024-27. Author ABK was also supported by an FPU predoc-toral grant and a mobility grant, both of the Ministry of Universities of the Spanish Government, and by the Hungarian grant of the National Research, Development and Innovation Office RRF-2.3.1-21-2022-00006. We also want to thank the reviewers for their useful comments.

## Conflicts of interest

The authors declare no conflicts of interest.
References

[1] F. Anstett-Collin, L. Denis-Vidal, and G. Millérioux. "A priori identifiability: An overview on definitions and approaches". In: *Annual Reviews in Control* 50 (2020), pp. 139-149.

[2] J. Arino and S. Portet. "A simple model for COVID-19". In: *Infectious Disease Modelling* 5 (2020), pp. 309-315.

[3] X. Rey Barreiro and A.F. Villaverde. "Benchmarking tools for a priori identifiability analysis". In: *Bioinformatics* 39.2 (2023), btad065.

[4] G. Bastin and D. Dochain. *On-line Estimation and Adaptive Control of Bioreactors*. Process Measurement and Control. Elsevier, 1990.

[5] G. Bellu et al. "DAISY: a new software tool to test global identifiability of biological and physiological systems". In: *Computer Methods and Programs in Biomedicine* 88.1 (2007), pp. 52-61.

[6] A. Bicchi et al. "On the observability of mobile vehicles localization". In: *Theory and Practice of Control and Systems*. World Scientific Pub Co Inc, 1999, pp. 142-147.

[7] M. Bocher. "The Theory of Linear Dependence". In: *Annals of Mathematics* 2.1/4 (1900), pp. 81-96.

[8] M. Bongarti et al. "Alternative SIAR models for infectious diseases and applications in the study of non-compliance". In: *Mathematical Models and Methods in Applied Sciences* 32.10 (2022), pp. 1987-2015.

[9] F. Brauer, C. Castillo-Chavez, and Z. Feng. "Models for Influenza". In: *Mathematical Models in Epidemiology* 69 (2019), pp. 311-350.

[10] F. Brauer and C. Castillo-Chávez. *Mathematical Models in Population Biology and Epidemiology*. Springer, 2012.

[11] F. Brauer et al. "Some models for epidemics of vector-transmitted diseases". In: *Infectious Disease Modelling* 1.1 (2016), pp. 79-87.

[12] A.P. Browning et al. "Identifiability analysis for stochastic differential equation models in systems biology". In: *Journal of the Royal Society - Interface* 17 (2020), p. 20200652.

[13] M.A. Capistrán, M.A. Moreles, and B. Lara. "Parameter estimation of some epidemic models. The case of recurrent epidemics caused by respiratory syncytial virus". In: *Bulletin of Mathematical Biology* 71.8 (2009), pp. 1890-1901.

[14] D. Champredon, J. Dushoff, and D.J.D. Earn. "Equivalence of the Erlang-distributed SEIR epidemic model and the renewal equation". In: *SIAM Journal on Applied Mathematics* 78.6 (2018), pp. 3258-3278.

[15] Y. Chen and J. Zhang. "Dynamic analysis and optimal control of influenza model with discrete age structure". In: *Mathematical Methods in the Applied Sciences* 47.6 (2024), pp. 4260-4282.

[16] G. Chowell. "Fitting dynamic models to epidemic outbreaks with quantified uncertainty: A primer for parameter uncertainty, identifiability, and forecasts". In: *Infectious Disease Modelling* 2.3 (2017), pp. 379-398.

[17] G. Chowell et al. "Structural identifiability analysis of epidemic models based on differential equations: a tutorial-based primer". In: *Journal of Mathematical Biology* 87.79 (2023).

[18] R. Csuma-Kovács and G. Röst. "Reproduction numbers for epidemiological models with realistic age structure (RAS models)". In: *Advances in Continuous and Discrete Models* 2025.163 (2021).

[19] N. Cunniffe et al. *Observability, Identifiability and Epidemiology – A primer*. Springer, 2024.
[20] V.T. Huynh D.C. Huong and H. Trinh. "Design of event-triggered interval functional observers for systems with input and output disturbances". In: *Mathematical Methods in the Applied Sciences* 44.18 (2021), pp. 13968-13978.

[21] J.F. David et al. "A novel approach to modelling the spatial spread of airborne diseases: an epidemic model with indirect transmission". In: *Mathematical Biosciences and Engineering* 17.4 (2020), pp. 3294-3328.

[22] A. Dénes and A.B. Gumel. "Modeling the impact of quarantine during an outbreak of Ebola virus disease". In: *Infectious Disease Modelling* 4 (2019), pp. 12-27.

[23] L. Denis-Vidal, G. Joly-Blanchard, and C. Noiret. "Some effective approaches to check the identifiability of uncontrolled nonlinear systems". In: *Mathematics and Computers in Simulation* 57 (2001), pp. 35-44.

[24] S. Diop and M. Fliess. "Nonlinear observability, identifiability, and persistent trajectories". In: *Proceeding of the 30th on Decision and Control*. 1981.

[25] R. Dong et al. "Differential Elimination for Dynamical Models via Projections with Applications to Structural Identifiability". In: *SIAM Journal on Applied Algebra and Geometry* 7.1 (2023).

[26] N.D. Evans et al. "The structural identifiability of the susceptible infected recovered model with seasonal forcing". In: *Mathematical Biosciences* 194 (2005), pp. 175-197.

[27] M. Fang and P. -A. Bliman. "Modelling, Analysis, Observability and Identifiability of Epidemic Dynamics with Reinfections". In: *2022 European Control Conference (ECC)* (2022), pp. 278-283. URL: https://api.semanticscholar.org/CorpusID:246275908.

[28] Z. Feng, C. Castillo-Chavez, and A.F. Capurro. "A Model for Tuberculosis with Exogenous Reinfection". In: *Theoretical Population Biology* 57.3 (2000), pp. 235-247.

[29] Z. Feng and J.X. Velasco-Hernández. "Competitive exclusion in a vector-host model for the dengue fever". In: *Journal of Mathematical Biology* 35 (1997), pp. 523-544.

[30] Z. Feng et al. "Computation of $\mathcal{R}$ in age-structured epidemiological models with maternal and temporary immunity". In: *Discrete and Continuous Dynamical Systems - B* 21.2 (2016), pp. 399-415.

[31] R. Gani and S. Leach. "Transmission potential of smallpox in contemporary populations". In: *Nature* 414 (2001), pp. 748-751.

[32] J.P. Gauthier, H. Hammouri, and S. Othman. "A Simple Observer for Nonlinear Systems – Applications to Bioreactors". In: *IEEE Transactions on Automatic Control* 37.6 (1992), pp. 875-880.

[33] D.K. Geller and I. Klein. "Angles-Only Navigation State Observability During Orbital Proximity Operations". In: *Journal of Guidance, Control, and Dynamics* 37.6 (2014), pp. 1976-1983.

[34] M. Gevers et al. "Identifiability and excitation of linearly parametrized rational systems". In: *Automatica* 63 (2016), pp. 38-46.

[35] S.T. Glad and L. Ljung. "On global identifiability for arbitrary model parametrizations". In: *Automatica* 30.2 (1994), pp. 265-276.

[36] S.T. Glad and L. Ljung. "Parametrization of nonlinear model structures as linear regressions". In: *IFAC Proceedings Volumes* 23.8/3 (1990), pp. 317-321.

[37] G.C. Goodwin and R.L. Payne. *Dynamic System Identification: Experiment Design and Data Analysis*. Academic Press, 1977.

[38] F. Hamelin et al. "About the identifiability and observability of the SIR epidemic model with quarantine". In: *IFAC-PapersOnLine* 56.2 (2023), pp. 4025-4030.

[39] R. Hermann and A.J. Krener. "Nonlinear Controllability and Observability". In: *IEEE Transactions on Automatic Control* 22.5 (1977), pp. 728-740.
[40] S. Holmsen, S. Eidnes, and S. Riemer-Sørensen. "Pseudo-Hamiltonian system identification". In: *Journal of Computational Dynamics* 11.1 (2024), pp. 59-91.

[41] B. Ivorra et al. "Mathematical modeling of the spread of the coronavirus disease 2019 (COVID-19) taking into account the undetected infections. The case of China". In: *Communications in Nonlinear Science and Numerical Simulations* 88 (2020), p. 105303.

[42] R. Jain, S. Narasimhan, and N.P. Bhatt. "A priori parameter identifiability in models with non-rational functions". In: *Automatica* 109 (2019), p. 108513.

[43] G. Jerónimo, M. Pérez Millán, and P. Solernó. "Identifiability from a Few Species for a Class of Biochemical Reaction Networks". In: *Bulletin of Mathematical Biology* 81 (2019), pp. 2133-2175.

[44] Q. Jiang et al. "A tuberculosis model with early and late latency, imperfect vaccination, and relapse: An application to China". In: *Mathematical Methods in the Applied Sciences* 46.9 (2023), pp. 10929-10946.

[45] R.E. Kalman. "A New Approach to Linear Filtering and Prediction Problems". In: *Journal of Basic Engineering* 82.1 (1960), pp. 35-45.

[46] R.E. Kalman. "On the General Theory of Control Systems". In: *IFAC Proceedings Volumes* 1.1 (1960), pp. 491-502.

[47] S. Kepley and T. Zhang. "A constructive proof of the Cauchy-Kovalevskaya theorem for ordinary differential equations". In: *Journal of Fixed Point Theory and Applications* 23.7 (2021), pp. 1-23.

[48] W.O. Kermack and A.G. McKendrick. "A contribution to the mathematical theory of epidemics". In: *Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences* 115.772 (1927), pp. 700-721.

[49] H.K. Khalil. *Nonlinear Systems (3rd Edition)*. Prentice Hall, 1996.

[50] 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.

[51] M. Kretzschmar, P.F.M. Teunis, and R.G. Pebody. "Incidence and reproduction numbers of pertussis: Estimates from serological and social contact data in five European countries". In: *PLoS Medicine* 7.6 (2010), e1000291.

[52] A.B. Kubik et al. "About identifiability and observability for a class of dynamical systems". Preprint in arXiv: https://doi.org/10.48550/arXiv.2408.11482. 2024.

[53] Y. Liu and J. Rocklöv. "The effective reproductive number of the Omicron variant of SARS-CoV-2 is several times relative to Delta". In: *Journal of Travel Medicine* 29.3 (2022), taac037.

[54] 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". In: *Journal of Travel Medicine* 28.7 (2021), taab124.

[55] Y.R. Liyanage et al. "A Tutorial on Structural Identifiability of Epidemic Models Using StructuralIden-tifiability.jl". Preprint in arXiv: https://doi.org/10.48550/arXiv.2505.10517. 2025.

[56] D.G. Luenberger. "Observers for Multivariable Systems". In: *IEEE Transactions on Automatic Control* 11.2 (1966), pp. 190-197.

[57] N. MacDonald. *Time Lags in Biological Models*. Springer Nature, 1978.

[58] A. Martinelli. *Observability: A New Theory Based on the Group of Invariance*. Society for Industrial and Applied Mathematics, 2020.

[59] G. Massonis, J.R. Banga, and A.F. Villaverde. "Structural identifiability and observability of compartmental models of the COVID-19 pandemic". In: *Annual Reviews in Control* 51 (2021), pp. 441-459.

[60] MathWorks. *Help Center - lsqcurvefit*. https://www.mathworks.com/help/optim/ug/lsqcurvefit.html.
[61] M. Medvedeva et al. "Direct estimation of SIR model parameters through second-order finite differences". In: *Mathematical Methods in the Applied Sciences* 44.5 (2021), pp. 3819-3826.

[62] A. Ovchinnikov, G. Pogudin, and P. Thompson. "Input-output equations and identifiability of linear ODE models". In: *IEEE Transactions on Automatic Control* 68.2 (2023), pp. 812-824.

[63] A. Ovchinnikov et al. "Computing all identifiable functions of parameters for ODE models". In: *Systems & Control Letters* 157 (2021), p. 105030.

[64] N. Parolini et al. "Modelling the COVID-19 epidemic and the vaccination campaign in Italy by the SUIHTER model". In: *Infectious Disease Modelling* 7.2 (2022), pp. 45-63.

[65] L. Plötzke et al. "Revisiting the Linear Chain Trick in epidemiological models: Implications of underlying assumptions for numerical solutions". In: *Mathematics and Computers in Simulation* 239 (2026), pp. 823-844.

[66] Á.M. Ramos et al. "A simple but complex enough $\theta$-SIR type model to be used with COVID-19 real data. Application to the case of Italy". In: *Physica D: Nonlinear Phenomena* 421 (2021), p. 132839.

[67] Á.M. Ramos et al. "Modeling the impact of SARS-CoV-2 variants and vaccines on the spread of COVID-19". In: *Communications in Nonlinear Science and Numerical Simulation* 102 (2021), p. 105937.

[68] A. Rapaport and D. Dochain. "Interval observers for biochemical processes with uncertain kinetics and inputs". In: *Mathematical Biosciences* 193.2 (2005), pp. 235-253.

[69] W. Rudin. *Principle of Mathematical Analysis (3rd Edition)*. McGraw Hill, 1976.

[70] M.P. Saccomani, S. Audoly, and L. D'Angiò. "Parameter identifiability of nonlinear systems: the role of initial conditions". In: *Automatica* 39 (2003), pp. 619-632.

[71] E. Scharbarg et al. "From the hospital scale to nationwide: observability and identification of models for the COVID-19 epidemic waves". In: *Annual Reviews in Control* 50 (2020), pp. 409-416.

[72] R. Seck et al. "An age structured mathematical model for studying Malaria transmission dynamics. Applications to some areas of Senegal". In: *Mathematics and Computers in Simulation* 229 (2025), pp. 392-408.

[73] R. Seck et al. "An optimal control model to design strategies for reducing the spread of the Ebola virus disease". In: *Mathematical Biosciences and Engineering* 19.2 (2022), pp. 1746-1774.

[74] E. Tunali and T.-J. Tarn. "New results for identifiability of nonlinear systems". In: *IEEE Transactions on Automatic Control* 32.2 (1987), pp. 146-154.

[75] A.F. Villaverde, A. Barreiro, and A. Papachristodoulou. "Structural Identifiability of Dynamic Systems Biology Models". In: *PLOS Computational Biology* 12.10 (2016), e1005153.

[76] E. Walter and L. Pronzato. *Identification of parametric models from experimental data*. Springer, 1997.

[77] L. Wang, R. Ortega, and A. Bobtsov. "Observability is sufficient for the design of globally exponentially stable state observers for state-affine nonlinear systems". In: *Automatica* 149 (2023), p. 110838.

[78] Z.S.Y. Wong et al. "A systematic review of early modelling studies of Ebola virus disease in West Africa". In: *Epidemiology & Infection* 145.6 (2017), pp. 1069-1094.

[79] World Health Organization. *Consensus document on the epidemiology of severe acute respiratory syndrome (SARS)*. https://iris.who.int/handle/10665/70863.

[80] World Health Organization. *Ebola disease caused by Sudan ebolavirus - Uganda*. https://www.who.int/emergencies/disease-outbreak-news/item/2022-DON423.

[81] World Health Organization. *WHO Coronavirus (COVID-19) Dashboard*. https://covid19.who.int/.