One of the main problems faced in the mathematical modeling of the coronavirus epidemic has been the lack of quality data. In particular, it is estimated that a large number of cases have been unreported, especially those of asymptomatic patients. This is mostly due to the strong demand of tests required by a relatively complete report of the infected cases. Although the countries/regions that have managed to control the epidemic have been precisely those that have been able to develop a great capacity of testing, this has not been achieved in most of the situations.
In general, the number of unreported patients has been estimated by extrapolating data from the reported cases assuming that both numbers vary proportionally. However, this view can be criticized (al least) in two directions:
-It does not consider the dynamical role of the unreported cases in the evolution of the epidemic and, by extension, in the number of reported cases; -it does not allow a possible variation of the proportion between the numbers of reported and unreported cases when the conditions of social distancing remain unchanged.
In this work, we address these points incorporating the unreported cases into the modeling. In a first introductory section, we discuss from a mathematical perspective what happens when the testing capacity is very low. Using a very simple argument we show that, in this context, the dynamics of the disease is actually governed by the growth of the number of unreported cases, and that of the reported patients becomes much smaller (in a progressive way) than that of the total cases. Although the discussion is presented in a context of extreme (hence "ideal") conditions, it certainly illustrates how essential is to consider the role of the unreported cases to analyze the global evolution of the epidemic.
Several mathematical models have been proposed to deal with the unreported cases and their role in the progression of the epidemic. In the second section of this work we address one of them, with the acronym SIRU, recently proposed by Zhihua Liu, Pierre Magal, Ousmane Seydi and Glen Webb [10, 11, 12] . After a brief presentation of the model and an explanation of why it is more adjusted to the current epidemic, we proceed to establish a series of structural results. Although this does not completely close the study of the qualitative properties of the underlying differential equations, we can already glimpse an analogy with those of the classical SIR model. However, we stress an important difference: in the SIRU model, the curves that appear do not necessarily have a single peak. Specifically, we exhibit a simple method to detect parameters that give rise to curves with at least two peaks.
The model of Liu, Magal, Seydi and Webb has already been used to describe the evolution of the epidemic in various countries (China, South Korea, the United Kingdom, Italy, France and Spain). In the last section of this work, we implement this modeling to the global Chilean scenario using the official COVID-19 data provided by the Chilean government [13] . However, unlike [10, 11, 12] , our implementation is novel in that it uses a variable transmission rate for the disease, which is more pertinent according to the local epidemiological evolution.
This work concludes with a section of general conclusions in which some future lines of research are also described.
1 On the dynamics with low testing capacity
We denote by U the number of positive cases that one would expect to detect within a certain unit of time in some region using the maximal capacity of testing that is available. 3 Suppose that the epidemic has evolved to a point where the number of detected cases in this unit time is systematically very close to this threshold U . The reported contagion curve is then taking a plateau shape, but a question immediately arises: has the contagion process entered into an stationary phase -with a reproduction rate R equal to 1 or slightly higher-, or is a significant number of positive cases being indetected?
It is impossible to answer a priori to this question. However, it is very unlikely that the process enters a stationary regime exactly when the threshold is reached. We argue below that, if the entry into an stationary regime has not occurred, then not only a fixed "proportion" of the number of cases is being undetected, but this is the case of a much larger amount; more precisely, the latter follows an exponential rate growth.
To simplify the discussion, our unit time will be equal to the total period of the disease. For an instant i, we denote by c i the number of reported positive cases, and by C i that of total cases. We assume that, while c i remains below the threshold, it assumes values very close to it in the future evolution. We write this as follows:
where v is relatively small compared to U .
We will also assume that we are not in a stationary regime (that is, the apparent stationarity is actually a consequence of a default of testing). Therefore, the value of C i is substantially greater than the threshold and, consequently, than c i . At the instant i, the number of cases that are not being detected is C i − c i . These individuals will have a higher social activity than those detected as infected (since the latter will be quarantined). Therefore, the average number of new infected individuals by these undetected agents will be a valueR strictly larger than R.
Thus, on average, the detected infected individuals of time i generate R new infections at time i + 1, while those undetected at time i (whose number is C i − c i ) each generateR new infections. Therefore, the following inequality holds:
Since c i ≥ U − v and c i+1 ≤ U , the above implies
By simple recurrence, for an initial time i 0 and all n ≥ 1, this gives
In other words,
It is natural to expect that the term v/(R − 1) is (very) small with respect to C i 0 − c i 0 . Indeed, on the one hand,R does not approximate indefinitely to 1 (since the undetected infected individuals do not significantly reduce their social activity); on the other hand, our analysis deals with times for which C i is very (although, a priori, not exponentially) greater than c i . (A more robust argument consists in choosing not only a single initial time i 0 , but to implement the previous inequality along a sequence of consecutive times and finally average along these inequalities.)
Assuming all of the above, the conclusion is clear: the number of infected people C i 0 +n is much higher than the number c i 0 +n of reported infected individuals. Indeed, the difference between the two is bounded from below by an exponential of ratioR, while that the evolution of the detected cases is governed by the rate R, which is strictly less thanR. So, the curve that we see (that of the values c i+n ) is not only very far from the real one (that of C i+n ), but the difference between them has an exponential growth that we are not perceiving.
The "toy argument" above shows something evident: if we are not able to follow the evolution of an epidemic through appropriate mass testing, then we lose track of the infectious curve. In more sophisticated terms, what it reveals is that as long we are not aware that the reproduction rate R is strictly smaller than 1, the dynamical system of the epidemic moves in a regime of either slight exponential growth or, at least, of instability. In such a regime, small variations of the initial conditions can lead to exponential explosion. Now, to a large extent, these initial conditions are provided by official data. However, if these move around the maximum of what the system can detect, then we can hardly know how accurate they are and, therefore, whether we really are in a stationary situation or whether we have advanced to an exponential explosion of cases that we are not perceiving. For this reason, for each instance in which the number of positive detected cases becomes nearly constant (that is, when the curve of reported cases begins to acquire a plateau shape), it seems reasonable to apply a substantive increase in the number of tests (both in quantity and spectrum). If this leads to a significant increase in the number of positive cases, then most likely this would mean that, actually, the regime was not stationary, but was simply exceeding a threshold above which a significant amount of infections cannot be reported. We will return to this point in the general conclusions of this work.
The key argument in the preceding section is that unreported patients have a greater dynamic role than those that are reported (since the former do not enter into quarantine), and therefore they contribute more importantly to the epidemic. However, in the traditional SIR model, both types of patients are part of the same compartment. In [10] , Liu, Magal, Seydi and Webb solve this problem by separating them into two compartments. Denoting respectively by S, I, R and U the susceptible individuals, infected individuals who do not yet have symptoms (and are at incubation stage), reported infected individuals, and unreported (either asymptomatic or low symptomatic) infected individuals, they consider the following diagram flux: The differential equations attached to the diagram above and that govern the dynamics of the epidemic are the following:
where t ≥ t 0 corresponds to time, with t 0 being the starting date for the study (as in [10, 11, 12] , in the implementation, we will consider the time t 0 corresponding to the beginning of the epidemic). Although this system of differential equations makes perfect sense when prescribing any initial condition, in epidemiological modeling one is naturally lead to use data of the following type:
The parameters used in the model are described in the Table below. In particular, note that ν = ν 1 + ν 2 . In addition, all the parameters that are considered τ, ν, ν 1 , ν 2 , η are positive.
Interpretation
Transmission rate of the disease. 1/ν Average time during which the infectious asymptomatic individuals remain in incubation. f
Fraction of asymptomatic infected individuals that become reported infected. ν 1 = f ν Rate at which asymptomatic infected cases become reported. ν 2 = (1−f )ν Rate at which asymptomatic infected become unreported infected individuals (asymptomatic or mildly symptomatic). 1/η Average time during which an infected individual presents symptoms. Note that in the first of the equations of (1), namely
the role of I and U in the spread of the infection is the same. This is the essential point of the model: it gives the same dynamic role to those who are infected and do not have symptoms as to those who are not reported, because the latter do not go into quarantine and, in fact, have a similar social activity. Certainly, the model can be refined in many directions; for example, one could give different though still dynamical roles to I and U by attaching to them different positive parameters τ 1 = τ 2 (this could be justified by that U includes individuals with low symptomaticity who can practice self-care). However, it is already worth to visualize some of the main properties of this model and to implement it in specific situations following this general format.
We next consider the system of equations (1) in more detail. We will assume an epidemic situation, which is summarized by the condition
We will also consider the initial conditions of Liu, Magal, Seydi and Webb:
for a certain χ 2 > 0. (The precise value of the parameter χ 2 will be given in the next section; here we just retain the fact that it is positive.) For simplicity, our starting time will be t 0 := 0.
Theorem 2.1. In an epidemic situation (2) and starting with the initial conditions (3), the following three properties are fulfilled:
(i) For all time t > 0, the values of S(t), U (t), I(t) and R(t) exist, are positive and strictly smaller than P := S 0 + I 0 + U 0 (the total population);
(ii) S(t) converges to a certain positive limit value S ∞ as t → ∞, while
Proof. We first note that, due to (2) and (3), we have
Therefore, there exists ε > 0 such that S(t), I(t), R(t) and U (t) (are defined in [0, ε) and) are strictly positive. The arguments that follow are inspired by an observation contained in the classical book of Vladimir Arnold [1] .
(i) Suppose S vanishes, and let T > 0 be the first time this occurs. Let C > 0 be such that
Integrating between 0 and s < T , this gives Finally, suppose that I vanishes, and let T be the first moment this occurs. Then, U (t) ≥ 0 for all t ∈ [0, T ], and therefore
However, by integration, this again gives a contradiction, namely I(T ) ≥ I 0 e −νT .
We next show that neither I nor U explode (that is, none of them tends to infinity along an increasing sequence of times tending to a finite time T ). Indeed, if anyone does it in time T then, from the above, P ≥ S(t) ≥ 0 on [0, T ), and U, I are positive on this interval. Therefore, on [0, T ),
which implies by integration that (I + U )(t) ≤ (I 0 + U 0 ) e τ P t for t ∈ [0, T ). Letting t go to T , this contradicts the explosion.
To see that R does not explode, we proceed again by contradiction: if this occurs at time T , then from R = ν 1 I − ηR ≤ ν 1 I we deduce R(t) ≤ R 0 e ν 1 C for all t ∈ [0, T ) and C := max t∈[0,T ] I(t).
Finally, S does not explode because it is decreasing and positive.
In conclusion, S, I, R, U are all positive and do not explode. To prove that they are bounded from above by P , we introduce the equation of the deletted (removed) individuals from the system:
We have a constant population P = S(t) + I(t) + R(t) + U (t) + D(t), and since D > 0, we have that D(t) > 0 for all t > 0. Now, since S, I, R, U are positive for t > 0, we conclude that each of them must be strictly smaller than P .
(ii) First we show that S converges towards a positive limit. Let S ∞ be the limit of S (which exists because S is decreasing). Denoting c := min{η, ν}, we have
Since S = −τ S [I + U ], this implies
Integrating between 0 and t > 0, this gives
Thus, for all t > 0, we have S(t) ≥ S 0 e − τ c P , and therefore, S ∞ ≥ S 0 e − τ c P > 0.
Next, we simultaneously prove that I and U converge to 0 (the convergence of R to 0 will be then a consequence of the convergence of R/U to ν 1 /ν 2 proved in (iii) below). To do this, we first note that, since
it follows that there is a sequence of times t n → ∞ such that (I + U )(t n ) → 0. Otherwise, there would exist c > 0 such that (I + U )(t) ≥ c for all t > 0, which implies S /S ≤ −τ c, and therefore S(t) ≤ S 0 e −cτ t . However, for t large enough, this contradicts the inequality
Now, since ε > 0, we may fix T such that (I + U )(T ) ≤ ε/2 and S(t) ≤ S ∞ + ε/2 for all t ≥ T . We claim that (I + U )(t) ≤ ε for all t ≥ T . (Since ε > 0 was arbitrarily chosen, this concludes the proof of the convergence of I + U towards 0.) To prove this, note that, since
as we wanted to show.
(iii) To prove that R/U converges to ν 1 /ν 2 , we first remark that
Therefore,
In other words, if R/U is smaller (resp. greater) than ν 1 /ν 2 at a point t, then it is increasing (resp. decreasing) around this point.
We first prove that R/U cannot be equal to ν 1 /ν 2 at any point. To do this, we note that R 0 /U 0 = 0 = ν 1 /ν 2 . We define ϕ := ν 1 ν 2 − R U . Equality (4) then becomes
If T were the first instant at which R/U = ν 1 /ν 2 , then T would be the first zero of ϕ. By (i), there exists C > 0 such that ν 2 I/U ≤ C on [0, T ]. Since R 0 = 0, one has ϕ(t) > 0 for all positive but small-enough t. Choosing such a t smaller than T we obtain, on [t, T ),
By integration, this gives ϕ(T ) ≥ ϕ(t) e t−T , which contradicts the fact that ϕ(T ) = 0.
To prove the convergence of R/U towards ν 1 /ν 2 , which is equivalent to that of ϕ towards 0, we will use the fact (proven below) that I/U is bounded from below by a positive constant c. Assuming this, we have ϕ ϕ ≤ −ν 2 c, and hence,
which converges to 0 as T → ∞.
To conclude, we must show that I/U does not approach zero. For this, we begin by noting that
Since S ≥ S ∞ , if I/U is very small, then the derivative (I/U ) becomes positive, and therefore I/U grows. As a consequence, I/U cannot arbitrarily approximate 0.
There are several remarks to the proof above.
Remark 2.2. The proof above applies to more general initial conditions than (3): it only requires that the values S 0 > 0, I 0 > 0, R 0 ≥ 0 and U 0 ≥ 0 are such that I (t 0 ), R (t 0 ) and U (t 0 ) are all positive. However, note that, at this level of generality, in statement (iii) above the convergence of R/U towards ν 1 /ν 2 can occur from above, and even the quotient R/U can remain constant and equal to ν 1 /ν 2 throughout the whole evolution.
In the classical SIR model, the fact that the population of susceptibles converges towards a positive limit is often presented as a consequence of the so-called final size relation [5] . In our context, S not only depends on I, but also on U . For this reason, it is hard to expect such a simple relation, and this partly justifies the use of robust estimates in the preceding proof.
Remark 2.4. The convergence of R and U towards 0 must hold at a speed lower than e −ηt . Indeed, from the last two equations of the system one obtains
which can be rewritten in the form ν 2 (Re ηt ) = ν 1 (U e ηt ) . By integration, this gives
Therefore,
Since R/U converges to ν 1 /ν 2 , the product U e ηt must tend to infinity, thus showing our claim for U . The claim for R immediately follows from this.
Remark 2.5. Another observation is related to the proportion R/U . Since this number grows throughout the epidemic, the dynamic variability that we mentioned at the beginning of this work holds. Moreover, the convergence of R/U towards ν 1 /ν 2 tells us what should be the values of the parameters in the implementation. We will return to this point for a specific case in the next section.
Remark 2.6. Finally, we must mention that one of the basic results on the SIR model has not been incorporated into the theorem above, namely, that of the uniqueness of the peak of the curve of infected cases. In fact, this uniqueness is no longer valid for the SIRU model. Below we draw an example of a double peak for the reported and unreported cases curves. Note that though this occurs under initial conditions different from (3), it holds in a context in which the theorem is still valid, according to Remark 2.2. Examples of this type can be easily obtained. For the curve R, one starts by imposing the conditions R (t p ) = 0 and R (t p ) > 0 for a time t p > 0 different from the starting time. These conditions can be translated into conditions only on the values of S, I, R and U at that time. Then the SIRU equations are implemented in both directions of time around t p . This moment will hence correspond to a local minimum of the curve, necessarily located between two peaks.
The same argument allows to build examples in which the curve of infected cases has two peaks. Naturally, this is due to the presence of U in the derivative of I, which is an absent element of the SIR model. In fact, in this one, the condition I = 0 necessarily implies I < 0, as a straightforward computation shows. This occurs only at the peak of the curve, which corresponds to the time when herd immunity is achieved.
Our goal now is to compare the data available on COVID-19 in Chile until May 14 (2020) [13] with what is predicted by the SIRU model, so that we can estimate the evolution of the number of unreported cases throughout the period. We point out that a first estimate of unreported cases in (some regions of) Chile using the data of March 2020 and the SIRU model was carried out by Mónica Candia and Gastón Vergara-Hermosilla in [6] .
To begin with, let us point out that in outbreaks of influenza disease, the parameters τ, ν, ν 1 , ν 2 , η, as well as the initial values S 0 , I 0 and U 0 , are generally unknown. However, it is possible to identify them from specific time data of reported symptomatic cases. The cumulative number of infectious cases reported at a time t, denoted by CR(t), is given by
This is data that is openly available. Likewise, the cumulative number of unreported cases at a time t is
We note that, up to constants, these values coincide respectively with
and U (t) + η Following the scheme used by Liu, Magal, Seydi and Webb, we will assume that, at an early stage of the disease, CR(t) has an (almost) exponential form:
For simplicity, we will assume that χ 3 = 1. The values of χ 1 and χ 2 will then be adjusted to the accumulated data of cases reported in the early phase of the epidemic 4 using a classical least squares method (after passing to logarithmic coordinates). According to the above (see [10] for details), for numerical simulations, the initial time for the beginning of the exponential growth phase is fitted at
Again, for simplicity, we will identify the initial value S 0 to that of the total population of Chile (since there is no prior immunity against the virus). Once the values of ν, η and f are set, the conditions at the beginning of the disease are naturally fitted as
We recall that 1/ν corresponds to the average time during which patients are asymptomatic infectious. As in [10, 11, 12] , we will let this parameter be equal to 7 (days), and the same value will be used for the average time in which patients are reported or unreported infectious:
Note that although the incubation period has been reported as being slightly smaller, the lag in the delivery of the results of the PCR tests in Chile justifies this choice.
To implement the SIRU system we still need to establish a good value for the parameter f (the fraction of symptomatic cases that are reported). Once this is fitted, we will have the values of
and we will just need to fit the value of τ .
Before continuing, it is worth pointing out that in [10, 11, 12] there is no major discussion on the criterion used to establish the value of f in the different scenarios. Actually, a value issued by the medical counterpart is assumed as valid. (For example, f = 0.8 is considered for China.) In our modeling, we will use the work of Baeza-Yates [4] and that of Castillo and Pastén [7] , who use the case fatality rate of the disease (with the correct correction according to its duration; see [2, 9] ) to give estimates for the right number of infected individuals. 5 Since Baeza-Yates' argument is simpler and is not included in an academic publication, we borrow it below in a language closer to that of the SIRU model. As we will see, it yields a method to adjusting the value of f that can be used in almost all contexts.
Since we know that R/U converges towards ν 1 /ν 2 , we assume for this calculation that R/U is simply equal to ν 1 /ν 2 . In addition, we will argue in discrete units of time (in days). Then we have
If d is the average time of illness to death and M (n) the number of deaths in day n, then the case fatality rate L corresponds to
.
Moreover, the reported case fatality rate is
In the Chilean context, deaths in the period studied occurred within 9.4 days after the disease was reported. Adjusting d = 9, the computation of L R is made feasible from the data available in [13] .
Finally, regarding the natural case fatality ratio L of the disease, it is deduced from international studies that, after adjusting it to the age distribution of the Chilean population, it should vary between 0.2% and 1%, with a very high tendency to be around 0.6%. In summary, this gives a parameter f varying between 0.1 and 0.5, with a high tendency to be close to 1 3.4 ∼ 0.3.
Given the heterogeneity of the safeguard measures taken by the Chilean government, instead of directly applying the SIRU model, it became more pertinent to us to consider a variable transmission rate τ (as a function of time). To analyze the variation of the values of τ (t), we observe how the percentage of the Chilean population subjected to confinement has been changing, which is illustrated below: We hence propose a function τ (t) of the form
where the I i 's correspond to successive time intervals. For the Chilean context, according to the graph of quarantines illustrated above (Figure 3 ), the chosen time intervals are those described in Table 2 Following [10, 11, 12] , the value of the transmission rate in I 1 is fitted to
Then, the parameters µ i are chosen in such a way that the reported cumulated cases in the numerical simulation align with the data of the reported cumulative number of infections at time t according to [13] . This is implemented with various values of f , following the guidelines of the preceding subsection (specifically, we work with f = 0.2, f = 0.3 and f = 0.4). A summary is described in the following Table: Parameter Using the data of the cumulated reported cases available in [13] , we can finally proceed to the simulations, which are shown below. They illustrate numerical estimates for the curves of CR(t), CU (t), R(t) and U (t). 
The epidemic outbreak of the new human coronavirus COVID-19 was first detected in Wuhan, China, in late 2019. In Chile, the first case was reported on March 3, 2020, in Talca (Maule region). Since then, modeling the epidemic in the country has faced the problem of the low availability of disaggregated data [3] .
The first part of this work was inspired by the situation experienced in Chile during the month of April, when the number of reported cases stabilized around 450 per day. According to most specialists, this was not accurately representing the genuine epidemiological evolution. Although it is difficult to argue that our reasoning from the first section fully applies to this situation (in particular, because the number of daily tests during the period was variable), the explosive increase of cases that subsequently occurred should call us to reflection on the point. Nevetheless, it is clear that this increase was also due in part to the relaxation of the protection measures, which is reflected not only in the absolute (and exponential) increase in the number of cases, but also in the proportion of positive cases with respect to tests (the latter despite of the significant increase in the number of tests [8] ).
The first part of the work naturally led us to model the dynamics of the epidemic incorporating a compartment for unreported cases so that we could deal with their active role in the evolution. To do this, we used the SIRU model recently introduced/implemented by Liu, Magal, Seydi and Webb in [10, 11, 12] . Since the qualitative theory of the underlying differential equations has not yet been treated, we developed some essenttial points of it in the second section, thus rigorously establishing fundamental structure theorems. However, we pointed out an important difference with the SIR model, namely, the possibility of multiple peaks for the curves of infected, reported and unreported patients. For the future, it would be desirable to complete the qualitative description of the solutions of the equations with respect to the parameters, with an emphasis on the phenomenon of multiplicity of peaks. For example, until now we do not know whether there can be more than two peaks and/or whether there is any restriction on the relative position between them. Without any doubt, this discussion is relevant for the implementation of disease containment policies, since it is directly related to the recognition of the moment when the curves begin to definitively descend.
In the last section of the work, we implemented an extension of the model to the case of Chile. To do this, we considered a variable transmission rate, which is more appropriate according to the local reality. Our rate is coupled to the official statistics provided by the government in [13] , information that also allowed us to make the parametric afjustement. The last ingredient to launch the simulations was to fit the value of the fraction f of infected individuals that were reported during the period of study. For this we used the work of Baeza-Yates [4] and that of Castillo and Pastén [7] . Naturally, our modeling is consistent with their estimates. In particular, in Section 2.3 of [7] , the authors estimate the number of real cases for April 28 6 as 43095, which is remarkably close to the number of total cases that can be inferred from the plot (B) in Figure 4 obtained by considering f = 0.3. We point out that the prediction of Castillo and Pastén is recovered with complete accuracy through the method used in this work for the value f = 0.31791. According to our modeling, on April 28, a percentage of 31.791% of the infections was reported to health agencies, and consequently 68.209% of the infected cases was not reported.
Although the estimates above were obtained a posteriori, their complementarity puts us in a good position to use them for modeling the future evolution of the epidemic. However, working in this direction requires great caution. In particular, it would be necessary to consider a variable parameter f , since both the quantity and the criteria of the tests have been modified during late May. Despite this, we strongly believe that the basis for pursuing the implementation of the SIRU model are fullfilled, and it would be very useful to advance in a more compartmentalized and georeferenced implementation of it.