**Proof of Theorem 2.**

(a) Following the method in [24], with reference to Equations (5)–(8), we have the rate of appearance of new infections vector F(*x*) and the rate of transfer of individuals vector V(*x*):

$$\mathcal{F} = \begin{pmatrix} 0\\ \beta SI\\ 0\\ 0 \end{pmatrix} \text{ and } \mathcal{V} = \begin{pmatrix} \beta SI + \mu S - \mu\\ \delta E + \mu E\\ -\delta E + \alpha I + \mu I + \varepsilon I\\ -\alpha I + \mu R - \varepsilon I \end{pmatrix}.$$

Note that there are only two sub-classes that involve infected persons (i.e., *E* and *I*), meaning we have *F* and *V* as 2 × 2 matrices. Here, we count a new infection as only occurring in *E* with the rate *βSI* and do not count the rate *δE* in *I* as a new infection, since it is only the transition from *E* to *I*. Next, from the two vectors, we obtain two matrices: *F* = - 0 *β* 0 0 and *V* = - *δ* + *μ* 0 −*δ α* + *μ* + *ε* . Consequently, *V*−<sup>1</sup> = 1 *<sup>δ</sup>*+*<sup>μ</sup>* 0 *δ* (*δ*+*μ*)(*α*+*μ*+*ε*) <sup>1</sup> *<sup>α</sup>*+*μ*+*<sup>ε</sup>* and *FV*−<sup>1</sup> = *βδ* (*δ*+*μ*)(*α*+*μ*+*ε*) *β α*+*μ*+*ε* 0 0 , which gives rise to the effective reproduction number <sup>R</sup>*<sup>ε</sup>* <sup>0</sup> <sup>=</sup> *<sup>ρ</sup>*(*FV*−1) = *βδ* (*α*+*μ*+*ε*)(*δ*+*μ*).


**Theorem 3.** *The SEIR model in Equations (5)–(8) always has a trivial equilibrium* (*S*<sup>0</sup> <sup>∗</sup>, *E*<sup>0</sup> ∗, *I*<sup>0</sup> ∗, *R*<sup>0</sup> ∗), *while the non-trivial equilibrium* (*Se* ∗, *Ee* ∗, *Ie* ∗, *Re* ∗) *exists only if the effective reproduction number is greater than 1*, *i.e.,* <sup>R</sup>*<sup>ε</sup>* <sup>0</sup> <sup>=</sup> *βδ* (*α*+*μ*+*ε*)(*δ*+*μ*) > 1.

**Proof of Theorem 3.** It is obvious as a consequence of Theorems 1 and 2. -

Up to this point, we concluded that the threshold we found earlier (i.e., <sup>T</sup> *<sup>ε</sup>* ) is actually equivalent to the "true" effective reproduction number, <sup>R</sup>*<sup>ε</sup>* <sup>0</sup> (Theorem 2(c)). Here, we could actually find another threshold that has the same threshold value as <sup>R</sup>*<sup>ε</sup>* <sup>0</sup>. Remember that in the derivation of the basic reproduction number, we noticed that there are only two sub-classes involving infected persons, i.e., *E* and *I*. Then, we have *F* and *V* as 2 × 2 matrices. Here, we count new infections only in *E* with the rate *βSI* and do not count the rate *δE* in *I* as a new infection, since it is only the transition from *E* to *I*. However, if we count new infections in *E* with the rate *βSI* and do count the rate *δE* in *I* as a new infection, then we have:

$$\begin{aligned} \mathcal{F} &= \begin{pmatrix} 0 \\ \beta SI \\ \delta E \\ 0 \end{pmatrix} \text{and } \mathcal{V} = \begin{pmatrix} \beta SI + \mu S - \mu \\ \delta E + \mu E \\ \alpha I + \mu I + \varepsilon I \\ -\alpha I + \mu R - \varepsilon I \end{pmatrix}, \\\text{Then we will have:} \\\ F &= \begin{pmatrix} 0 & \beta \\ \delta & 0 \end{pmatrix} \\\ \text{and:} \\\ V &= \begin{pmatrix} \delta + \mu & 0 \\ 0 & \alpha + \mu + \varepsilon \end{pmatrix} \\\ \text{giving:} \\\ V^{-1} &= \begin{pmatrix} \frac{1}{\delta + \mu} & 0 \\ 0 & \frac{1}{\alpha + \mu + \varepsilon} \end{pmatrix}. \end{aligned}$$

Hence, we obtain:

$$\begin{aligned} FV^{-1} &= \begin{pmatrix} 0 & \frac{\beta}{\alpha + \mu + \varepsilon} \\ \frac{\delta}{\delta + \mu} & 0 \end{pmatrix}, \\ \text{which implies that:}\\ \rho(FV^{-1}) &= \sqrt{\frac{\beta\delta}{(\alpha + \mu + \varepsilon)(\delta + \mu)}}. \end{aligned}$$

Note that the last expression is actually the square root of the reproduction number, so that by referring to [24], we have an alternative threshold, <sup>A</sup>*<sup>ε</sup>* <sup>=</sup> *<sup>ρ</sup>*(*FV*−1) = R*<sup>ε</sup>* <sup>0</sup>, which is not a reproduction number but has the same threshold value, i.e., 1.

**Theorem 4.** *The non-endemic equilibrium* (*S*<sup>0</sup> <sup>∗</sup>, *E*<sup>0</sup> ∗, *I*<sup>0</sup> ∗, *R*<sup>0</sup> ∗) *of Equations (5)–(8) is asymptotically stable whenever* <sup>R</sup>*<sup>ε</sup>* <sup>0</sup> <sup>=</sup> *βδ* (*α*+*μ*+*ε*)(*δ*+*μ*) < 1 *and unstable otherwise.*

**Proof of Theorem 4.** Let us consider the non-endemic equilibrium (*S*<sup>0</sup> <sup>∗</sup>, *E*<sup>0</sup> ∗, *I*<sup>0</sup> ∗, *R*<sup>0</sup> ∗) = (1, 0, 0, 0). The Jacobian matrix at this point is given by: ⎡ ⎢ ⎢ ⎣ −*μ* 0 −*β* 0 0 −*δ* − *μ β* 0 0 *δ* −*α* − *μ* −*ε* 0 0 *α ε* − *μ* ⎤ ⎥ ⎥ ⎦, which has the polynomial characteristics: *a*4*λ*<sup>4</sup> + *a*3*λ*<sup>3</sup> + *a*2*λ*<sup>2</sup> + *a*1*λ* + *a*<sup>0</sup> = 0 with: *a*<sup>4</sup> = 1; *a*<sup>3</sup> = (4*μ* + *δ* + *α* − *ε*); *<sup>a</sup>*<sup>2</sup> <sup>=</sup> *<sup>μ</sup>*(*<sup>δ</sup>* <sup>+</sup> <sup>3</sup>*<sup>μ</sup>* <sup>+</sup> *<sup>α</sup>* <sup>−</sup> *<sup>ε</sup>*) + *δα* <sup>+</sup> <sup>2</sup>*μδ* <sup>−</sup> *δε* <sup>+</sup> <sup>2</sup>*μα* <sup>+</sup> <sup>3</sup>*μ*<sup>2</sup> <sup>−</sup> <sup>2</sup>*με* <sup>−</sup> *βδ*); *<sup>a</sup>*<sup>1</sup> <sup>=</sup> *<sup>μ</sup>*(*δα* <sup>+</sup> <sup>2</sup>*μδ* <sup>−</sup> *δε* <sup>+</sup> <sup>2</sup>*μα* <sup>+</sup> <sup>3</sup>*μ*<sup>2</sup> <sup>−</sup> <sup>2</sup>*με* <sup>−</sup> *βδ*) +(*μδα* <sup>−</sup> *μδε* <sup>+</sup> *<sup>μ</sup>*2*<sup>δ</sup>* <sup>+</sup> *<sup>μ</sup>*2*<sup>α</sup>* <sup>−</sup> *<sup>μ</sup>*2*<sup>ε</sup>* <sup>+</sup> *<sup>μ</sup>*<sup>3</sup> <sup>+</sup> *βεδ* <sup>−</sup> *βμδ*);

Clearly, *<sup>a</sup>*<sup>4</sup> <sup>&</sup>gt; 0 and *<sup>a</sup>*<sup>3</sup> <sup>&</sup>gt; 0. Furthermore, we have *<sup>a</sup>*<sup>0</sup> <sup>&</sup>gt; 0, provided <sup>R</sup>*<sup>ε</sup>* <sup>0</sup> < 1. The detail of the proof is as follows.

Proof of *a*<sup>0</sup> > 0

We need *<sup>μ</sup>*(*μδα* <sup>+</sup> *μδε* <sup>+</sup> *<sup>μ</sup>*2*<sup>δ</sup>* <sup>+</sup> *<sup>μ</sup>*2*<sup>α</sup>* <sup>+</sup> *<sup>μ</sup>*2*<sup>ε</sup>* <sup>+</sup> *<sup>μ</sup>*<sup>3</sup> <sup>+</sup> *βεδ* <sup>−</sup> *βμδ*) <sup>&</sup>gt; 0, which can be written as follows:

*<sup>μ</sup>*2(*δε* <sup>+</sup> *μδ* <sup>+</sup> *με* <sup>+</sup> *<sup>μ</sup>*2)+(*μδα* <sup>+</sup> *<sup>μ</sup>*2*<sup>α</sup>* <sup>+</sup> *βεδ* <sup>−</sup> *βμδ*)*<sup>μ</sup>* <sup>&</sup>gt; <sup>0</sup> *<sup>μ</sup>*(*δε* <sup>+</sup> *μδ* <sup>+</sup> *με* <sup>+</sup> *<sup>μ</sup>*2)+(*μδα* <sup>+</sup> *<sup>μ</sup>*2*<sup>α</sup>* <sup>+</sup> *βεδ* <sup>−</sup> *βμδ*) <sup>&</sup>gt; <sup>0</sup> *<sup>μ</sup>*(*δε* <sup>+</sup> *μδ* <sup>+</sup> *με* <sup>+</sup> *<sup>μ</sup>*<sup>2</sup> <sup>+</sup> *δα* <sup>+</sup> *μα*) <sup>&</sup>gt; *βδ*(*<sup>μ</sup>* <sup>−</sup> *<sup>ε</sup>*) *βδ*(*μ*−*ε*) *<sup>μ</sup>*(*ε*+*μ*+*α*)(*δ*+*μ*) < 1 *βδ* (*ε*+*μ*+*α*)(*δ*+*μ*) <sup>&</sup>lt; *<sup>μ</sup>* (*μ*−*ε*).

*<sup>a</sup>*<sup>0</sup> <sup>=</sup> *<sup>μ</sup>*(*μδα* <sup>−</sup> *μδε* <sup>+</sup> *<sup>μ</sup>*2*<sup>δ</sup>* <sup>+</sup> *<sup>μ</sup>*2*<sup>α</sup>* <sup>−</sup> *<sup>μ</sup>*2*<sup>ε</sup>* <sup>+</sup> *<sup>μ</sup>*<sup>3</sup> <sup>+</sup> *βεδ* <sup>−</sup> *βμδ*).

Here, we need *μ* > *ε* to make the inequality consistent, since all of the parameters are non-negative. Note that in this case, *<sup>μ</sup>* (*μ*−*ε*) <sup>&</sup>gt; 1; hence, if <sup>R</sup>*<sup>ε</sup>* <sup>0</sup> < 1, then the inequality *βδ*

(*ε*+*μ*+*α*)(*δ*+*μ*) <sup>&</sup>lt; *<sup>μ</sup>* (*μ*−*ε*) holds. Proof of *a*<sup>1</sup> > 0

Note that *a*<sup>1</sup> can be written in the form of *a*<sup>1</sup> = *a*<sup>11</sup> + *a*<sup>0</sup> with *a*<sup>11</sup> = *μ*(*δα* + 2*μδ* + *δε* + <sup>2</sup>*μα* <sup>+</sup> <sup>3</sup>*μ*<sup>2</sup> <sup>+</sup> <sup>2</sup>*με* <sup>−</sup> *βδ*). We also note that *<sup>a</sup>*<sup>11</sup> <sup>−</sup> *<sup>a</sup>*<sup>0</sup> <sup>=</sup> <sup>2</sup>*μ*<sup>3</sup> <sup>+</sup> *δμ*<sup>2</sup> <sup>+</sup> *<sup>μ</sup>*2*<sup>ε</sup>* <sup>+</sup> *<sup>μ</sup>*2*<sup>α</sup>* is positive, so that *a*<sup>11</sup> > 0, since we proved earlier that *a*<sup>0</sup> > 0. Furthermore, since both *a*<sup>11</sup> > 0 and *a*<sup>0</sup> > 0, then, consequently, *a*<sup>1</sup> > 0.

Proof of *a*<sup>2</sup> > 0

Note that *a*<sup>2</sup> = *a*<sup>21</sup> + *a*11/*μ* with *a*<sup>21</sup> = *μ*(*δ* + 3*μ* + *α* + *ε*) is clearly positive. Since *a*<sup>11</sup> > 0, then *a*<sup>2</sup> > 0. Therefore, we proved that all of the coefficients of the polynomial characteristics are positive. Consequently, from the Descartes rule of signs, all roots have negative real parts. This proves the stability of the disease-free equilibrium whenever R*ε* <sup>0</sup> < 1. -

**Theorem 5.** *If the endemic equilibrium* (*Se* ∗, *Ee* ∗, *Ie* ∗, *Re* ∗) *of Equations (5)–(8) exists (i.e., whenever* <sup>R</sup>*<sup>ε</sup>* <sup>0</sup> <sup>=</sup> *βδ* (*α*+*μ*+*ε*)(*δ*+*μ*) > 1*), then it is asymptotically stable*.

**Proof of Theorem 5.** The proof is analogous as before. -

Assuming that in the absence of CPT, the system has a large basic reproduction number (otherwise, the administration of CPT will not make sense), then when CPT is administered, we could compute the ratio of the effective reproduction number to the basic reproduction number as:

$$\mathcal{R}\_0^\varepsilon : \mathcal{R}\_0 = \frac{\not\delta}{(\iota + \mu + \iota)(\delta + \mu)} : \frac{\not\delta}{(\iota + \mu)(\delta + \mu)} = \frac{(\iota + \mu)}{(\iota + \mu + \iota)} < 1$$

Hence, clearly, <sup>R</sup>*<sup>ε</sup>* <sup>0</sup> : R0. Consequently, we have the following theorem.

\*\*Theorem 6.\*\* For the SEIR model in Equations (5)-(8), we have:

$$\text{S}\_{\varepsilon}^{\*} = \frac{1}{\mathcal{R}\_{0}^{\varepsilon}} > \frac{1}{\mathcal{R}\_{0}^{\*}} = \text{S}\_{0}^{\*} \text{ and } I\_{\varepsilon}^{\*} = \left(\mathcal{R}\_{0}^{\varepsilon} - 1\right)\frac{\mu}{\mathcal{P}} < \left(\mathcal{R}\_{0} - 1\right)\frac{\mu}{\mathcal{P}} = I\_{0}^{\*} \text{ with the difference}$$

$$\text{S}\_{\varepsilon}^{\*} - \text{S}\_{0}^{\*} = \frac{\varepsilon(\delta + \mu)}{\delta\beta} \text{ and } I\_{0}^{\*} - I\_{\varepsilon}^{\*} = \frac{\delta\mu\varepsilon}{(\delta + \mu)(\alpha + \mu)(\mu + \alpha + \varepsilon)}.$$

**Proof of Theorem 6.** The proof is obvious. -

Analysis for the case of the rate of CPT proportional to the number of recovered class can be conducted analogously. We do not present the results explicitly, since all proofs are similar to the one presented here. In the next section, we carried out several simulations to assess the impact of CPT both to the transient solution and to the equilibrium solution. The simulation was conducted by implementing the Runge–Kutta numerical scheme to determine the numerical solution of the system, and the results are presented numerically.

#### *3.2. Numerical Examples*

In this section, we present numerical examples to show the behavior of the SEIR model with and without the presence of convalescent plasma transfusion. The results, in general, support the analysis of the SEIR equilibrium solutions presented in the earlier section. For the numerical examples, we used the CPT parameters in Table 2 and the other parameters written in the respective resulting figures. The results are summarized in the figures that follow.

**Table 2.** Parameter values used for different scenarios of CPT implementation.


Figure 2 shows a graph of the disease dynamics predicted by the SEIR model for specific parameters (written in the figure caption) with a high basic reproduction number (approximately 2.3). It suggests that, eventually, the disease will endemic to a certain level with a stable number of infectives (approximately 2% of the total population). Figure 3 shows the changes in the graph when CPT is utilized to cure patients. It shows that, for a relatively high rate of CPT intervention, it can reduce the basic reproduction number to the effective reproduction number as low as 0.97, which leads to a stable disease-free equilibrium.

**Figure 2.** A graph of the SEIR model without CPT showing the susceptible and recovered classes (**a**) and the exposed and infective classes (**b**). The SEIR parameters were *μ* = 1/75, *β* = 0.95, *δ* = 0.90, and α = 0.40, with initial values of *S*<sup>0</sup> = 0.99, *E*<sup>0</sup> = 0, *I*<sup>0</sup> = 0.01, and *R*<sup>0</sup> = 0. The resulting equilibrium was approximately *S* = 44%, *E* = 1%, *I* = 2%, and *R* = 53% with the basic reproduction number R<sup>0</sup> = 2.265.

**Figure 3.** A graph of the SEIR model with CPT, showing the susceptible and recovered classes (**a**) and the exposed and infective classes (**b**). The SEIR parameters were *μ* = 1/75, *β* = 0.95, *δ* = 0.90, and α = 0.40, with initial values of *S*<sup>0</sup> = 0.99, *E*<sup>0</sup> = 0, *I*<sup>0</sup> = 0.01, and *R*<sup>0</sup> = 0. The CPT rate was assumed to be ε = 0.55. The resulting equilibrium was approximately *S* = 99%, *R* = 0.05%, and the remaining *E* and *I* were nearly zero. In this case, the resulting effective reproduction number was R<sup>0</sup> = 0.97 (less than 1).

Figure 3 also suggests that early application of CPT significantly reduces the risk of disease outbreak. In the early transmission of COVID-19, this was clearly not the case, since COVID-19 is a new and emerging disease; hence, the availability of CP was almost null in the beginning. The figure also suggests a practical consequence of creating a convalescent plasma bank. Now, looking closer at the graph in Figure 2, during the first 200 time steps, we have the graph in Figure 4. It can be seen that there were already many recovered patients; hence, the availability of CP stock may be justified. Suppose that at time *t* = 200, the health authority begins to apply CPT as a curative method, then we have Figure 5, which shows that CPT significantly drove the disease cases down to zero, eventually. This is among the promising findings suggested by the SEIR model.

**Figure 4.** A graph of the SEIR model without CPT (top figures) showing the susceptible and recovered classes (**a**) and the exposed and infective classes (**b**) as in Figure 2 but with a shorter time horizon. (**c**,**d**) With CPT. The bottom figures show similar graphs for the SEIR model with CPT as in Figure 3 but with a shorter time horizon.

Figure 5 shows a graph of the SEIR model as in Figure 2, assuming that in the beginning (i.e., during the time interval (0,200)), the health authority takes the "do nothing" decision in controlling the disease (blue and red circles), and then, at time *t* = 200, begins to implement CPT with a relatively high rate of implementation (approximately half of the infectives are given CPT, proportional to the number of infectives with ε = 0.55). The black dots reveal that the intervention quickly pulls the number of infectives to zero (b) while, at the same time, pushes the number of susceptibles upward (a).

**Figure 5.** Graphs of the SEIR model assuming the "do nothing" decision during the time interval (0,200), followed by the implementation of CPT with the rate proportional to the number of infectives. The graph in (**a**) shows the dynamics of the susceptibles, while the graph in (**b**) shows the dynamics of the infectives.

Figure 6 shows the effects of different rates of CPT on decreasing the number of infectives (hence, the height of the infective peak). Here, we assumed a scenario in which CPT is conducted with the rate proportional to the infectives (implicitly assuming an abundance of CP bloods). Figure 7 shows the effect of different scenarios (see Table 1) on the decrease in the infective numbers over time. All of the figures assumed that there was no limit for the health authority to set CPT rates, except in Figure 7f, in which it was assumed that the maximum CPT rate was the Maxserv (response function number 6 in Table 2).

**Figure 6.** Graphs showing the effects of different rates of CPT on decreasing the number of infectives (i.e., lowering the height of the infective peaks) (**a**,**b**) and, consequently, increasing the number of remaining susceptibles (**c**,**d**). The other parameters were the same as in Figure 2.

In Figure 7f, if the Maxserv is unbounded, then the result is the same as in Figure 3b, which is unlikely in reality. Figure 8 shows examples of the Maxserv graphs used in Figure 7. The effect of various Maxserv response functions on the numbers of infectives and the exposed population is shown in Figure 9.

In all of the numerical simulations above, we assumed that the effect of CPT was on increasing patient survival. The recent findings reported in [26] support the assumption we used in this paper. They carried out a meta-analytical approach to collect and analyze the daily survival data from all controlled studies that reported Kaplan–Meier survival plots. The authors showed that CPT contributes to improving the symptomatology and viral clearance. Furthermore, they pointed out that the aggregate Kaplan–Meier survival plot in their study revealed a good agreement pattern among all different studies in which CPT was generally associated with greater patient survival [26].

The historical evidence shows the promising results of applying the therapeutic treatment of CPT for critical patients infected by contagious diseases such as COVID-19. A more recent study provided strong evidence that if the convalescent plasma is transfused into patients within three days of the onset of illness, a 41% lower risk of death compared to patients transfused four or more days after onset of illness is demonstrated [26,27]. This remarkable result highlights an important role for the timely use of convalescent plasma

transfusion. We discussed all such situations in terms of mathematical modeling with real data collected from a secondary source in the numerical analysis.

**Figure 7.** Plots of the effect of various CPT scenarios on the dynamics of the infective (red) and exposed (green) classes, with the assumption that the CPT intervention was carried out from the beginning of the pandemic. The scenarios and parameters were the same as in Table 2. The response functions and the parameters in the response functions for the graphs in (**a**–**f**) are presented in Table 2. See also Figure 2 as the reference for the other epidemiological parameters. Note that if Maxserv is unbounded, then the result is the same as Figure 3b.

**Figure 8.** The graph in (**a**) shows a plot the of mass action term of *εI*(*t*)*R*(*t*) in Figure 7c. The graph in (**b**) shows a plot of the numerical response min(*εI*(*t*)*R*(*t*), Maxserv) used in Figure 7f with ε = 1 and Maxserv = 0.01.

**Figure 9.** Plots of the effect of various Maxserv CPT scenarios on the dynamics of infective (red) and exposed (green) classes, with the assumption that the CPT intervention is conducted from the beginning of the pandemic. The scenarios and parameters for the graphs in (**a**–**d**) are as in Table 2.

In this research, a standard analytical evaluation, from proving the existence of equilibria, their local stability, and their relationships with the reproduction numbers, was only carried out for the model with the response function for the CPT rate proportional to the number of infected class. The numerical simulations showed the consistency of other forms of the response functions with the findings in real phenomena such as reported in [26,27]. However, to obtain more prudent results, it is necessary to undertake a complete mathematical analysis of all proposed response functions that will presumably provide more in-depth insight and better implementation. Well-posedness of solutions and the threshold criteria for the global stability of equilibria should certainly be sought [28]. On the contrary, the numerical solutions presented in this paper were obtained by the RK45, although according to the current findings [29,30], there is a better scheme to produce a dynamically consistent numerical solution. Other future research could focus on the application of optimal control theory for disease prevention and control based on different nonlinear response functions for the therapeutic rate of CPT. In particular, a study concentrating on applying a wider class of control variables to formulate an optimal control problem in order to find better control and preventive strategies in CPT implementation would be beneficial as found in [31].

#### **4. Conclusions**

We presented a continuous SEIR epidemic model considering the effect of an intervention using convalescent plasma transfusion (CPT) to the infected class. We analyzed the model using the standard procedure and found the trivial and non-trivial equilibrium points of the system including their stability and relation to the basic reproduction number. In general, the effect of the application of CPT on the individual level resulted in a shorter time of infection and a higher survival rate for infected individuals that received CPT. Furthermore, we showed that at the population level, it could also decrease the peak of the outbreak as well as the length of the epidemic period. In this case, the decrease in the infection peak indicated the good effect of the use of CPT, which may eventually decrease

the burden of COVID-19 transmission. The model presented here is still simple in terms of biological and epidemiological complexity; hence, further refinement of the model is still needed to obtain a more realistic model and a more accurate prediction.

In this paper, we proposed various functional forms/numerical responses that could be used to model the CPT rate, but not all were evaluated analytically. These numerical responses can be considered as a parametric viewpoint of control strategy. Hence, analytical investigation regarding the use of these various functional forms is also worthy to explore the robustness of the results presented here and to generate some possible epidemiological precautions not explored in the current paper. The model in this paper also only assumed a single strain or variant. In reality, viruses are always changing and mutating, and this could cause a new strain or variant. How would the results here be affected by such phenomena? Another question important for future research involves finding optimal CPT strategies that minimize both the number of infections and the related costs of CPT implementation. Our future research will focus on the application of optimal control theory for disease prevention and control based on different nonlinear response functions for the therapeutic rate of CPT. In particular, we will concentrate on applying a wider class of control variables to formulate an optimal control problem in order to find better control and preventive strategies for implementation of CPT.

**Author Contributions:** Conceptualization, H.H. and A.K.S.; methodology, H.H., A.K.S. and M.H.A.B.; software, A.K.S.; validation, H.H. and R.R.; formal analysis, A.K.S.; investigation, R.R.; resources, R.R.; data curation, R.R.; writing—original draft preparation, H.H.; writing—review and editing, H.H. and M.H.A.B.; visualization, A.K.S.; supervision, A.K.S.; project administration, H.H.; funding acquisition, H.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Indonesian Government through the scheme "Penelitian Kompetitif Nasional dan Penelitian Dasar", contract number 1207/UN6.3.1/PT.00/2021.

**Institutional Review Board Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.
