**Abstract**

This chapter proposes a theoretical reflection on the ethical deliberation process in the allocation of beds and respirators, in the light of the Theory of Health Justice, the Accountability for Reasonableness approach and the principle of health equity of the Brazilian Unified Health System (SUS, as per its Portuguese acronym), before the COVID-19 pandemic in Brazil. The pandemic has become a serious threat to health systems, as installed capacity has been exceeded whether in terms of material resources, equipment, technology and human resources. Thus, according to the theory of Accountability for Reasonableness, a fair and deliberative process aims to ensure resource allocation through limits and constraints (reasonableness), but government responsibility derived from human rights must be considered, allowing for health programming (accountability). Faced with this scenario, where the situation will often require us to make choices, this chapter intends to discuss the assumptions for ethical deliberation, taking into account the context of the act and its foreseeable consequences.

**Keywords:** Equity in Access to Health Services, Resource Allocation, Social Justice, Coronavirus Infections, Pandemic

#### **1. Introduction**

Throughout December 2019, in Wuhan city, Hubei province, Chinese health authorities identified pneumonia cases of unknown aetiology, which were subsequently attributed to the emergence of a new coronavirus (SARS-Cov-2) [1]. The World Health Organisation (WHO) has declared a state of public health emergency of international concern. Originally called 2019 - n - CoV, the infection caused by the new coronavirus was officially named COVID-19 on 11th March 2020, when it was declared a pandemic [2].

According to the Ministry of Health, in Brazil, until 28th April 2021, 14,521,289 confirmed cases of infection by COVID-19 had been notified, with 79,726 deaths [3]. This advance has led government agencies to adopt strategies to control the

spread of the virus, such as closing schools and shops, banning meetings, instituting teleworking, isolating people in their homes, and also banning international and even national travels [2].

The Brazilian Unified Health System (SUS, as per its Portuguese acronym) has always been at the front line of public health emergencies in Brazil and it has not been different in the context of facing the COVID-19 pandemic. It is one of the largest health systems in the world and the fight against the pandemic has integrated all levels of healthcare, based on its guidelines and principles described in the organic laws of health n° 8080 and n° 8142, dated 1990 [4, 5]. SUS must respond in a planned, systematised and equitable way to this new challenge imposed.

In an analysis of SUS, in its 30 years of existence, the system's lack of priority and underfunding, associated with the country's economic crisis and especially, the Constitutional Amendment 95 (EC-95/2016), which freezes the public budget for twenty years, were highlighted as factors detrimental to the system [6]. Thus, the system is overloaded and, faced with the new pandemic context, needs investment in material resources and equipment, as well as physical and human resources [7].

In the context of the pandemic in China, substantial regional disparities in the availability and accessibility of health resources were observed that could partially explain the low mortality rates (despite the large number of cases) in the more developed provinces, such as Zhejiang (zero deaths among 1,171 confirmed cases) and Guangdong (four deaths among 1,322 cases). Following analysis of this data, the Chinese government moved to rapidly build new local medical facilities in the less developed provinces, a key measure that helped control the epidemic, protect front-line health professionals and mitigate the severity of the disease in patients [8].

A reflection was introduced by health professionals from the United States on the need for resources to care for a large number of patients with COVID-19 after the declaration of the pandemic in the country, questioning how prepared the US hospitals were in relation to material and human resources in the face of uncontrolled transmission of the virus and consequent challenges in clinical decision making. In moments of crisis, health resources must be allocated in an ethical, rational and structured way, so as to benefit a greater number of patients. It is pointed out that the main conflict in decision making is the definition of ethical criteria for the allocation of mechanical ventilation equipment and beds in Intensive Care Units (ICU) [1].

We need to think about what ethical duties should be observed in this scenario, when the health system is being demanded beyond its capacity. How to allocate resources equitably, ensuring resources to all levels of care and regions? [2].

In pandemic contexts, there is a direct concern with equity in the use of resources as a moral and legal imperative to respect, protect and fulfil human rights in relation to health service delivery, based on international norms and standards. Nevertheless, in relation to the human rights-based approach, "even if attention is paid to the philosophical foundations linked to relevant rights recognised in international law, this does not provide guidance on how these relevant rights help determine priorities among claimants" [9].

General theories of justice alone do not provide adequate guidance for the problem of priority setting. This is because these theories are too broad to provide precise answers to the key distributional questions necessary for priority setting. However, the Accountability for Reasonableness approach, suggested by Daniels, guides prioritisation and is consistent with human rights due to its concern with justice and the social determinants in its process, which enables it to be applied to health systems [9].

*Ethical Deliberation in the Allocation of Respirators and Beds during the Covid-19 Pandemic… DOI: http://dx.doi.org/10.5772/intechopen.98589*

Accordingly, this chapter proposes a theoretical reflection on the process of ethical deliberation in the allocation of ICU beds and respirators, in the light of the Theory of Health Justice [10] and Accountability for Reasonableness approach [9], both proposed by the American philosopher and politician Norman Daniels; and the principle of SUS equity, before the COVID-19 pandemic in Brazil.


**5. Conclusions**

In this study, a new method is introduced which allows the derivation of DHR and BHR values from tower albedometer measurements. This method was applied to derive DHRs and BHRs over 20 tower sites, including both homogeneous and heterogenous land surfaces from the FLUXNET, SURFRAD and BSRN networks between the years 2012 and 2016. The retrieved DHR and BHR values were directly compared with the satellite albedo values, including CGLS, MODIS and MISR retrievals. The MODIS 16-day albedo products show the best agreement with in situ retrievals, whilst the MISR near-simultaneous measurements show a similar good agreement with in situ retrievals for a smaller time-window. The CGLS 30-day products have larger biases than MODIS and MISR products. Overall, the direct intercomparison with tower albedometer derived values shows a better match over the homogeneous sites than the heterogenous sites. The agreement between tower and satellite retrieved DHR values are better than the BHR values. This is because DHRs are only measured at local solar noon, whereas BHRs are derived from measurements at all possible solar zenith angles.

A surface albedo upscaling method, for tower FoV albedos to coarse resolutions, is described. This method employs atmospherically corrected BRFs from high-resolution EO alongside coarseresolution albedos predicted from a MODIS BRDF climatology over a larger area as inputs. The high-resolution albedo values are retrieved from the MODIS BRDF derived albedo-to-nadirreflectance ratios. This method was applied to upscale tower measured DHR and BHR values to 1-km resolutions and compared with the CGLS products. These results imply that this surface albedo upscaling strategy can be applied to both homogeneous and heterogenous surfaces if the optimal sample size for optimising this upscaling is known. For example, one of the sites, where 3 × 3 km were compared, shows that the pixel to the north-west of where the tower is located yields a better correlation than the pixel containing the tower, which appears to be associated with the observation that the land cover changes in the south-east. This will be explored in more detail in future.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2072-4292/11/6/644/s1, Figure S1: CGLS, MODIS and MISR DHR products compared with tower derived DHRs. Figure S2: CGLS, MODIS and MISR BHR compared with tower derived BHRs. Figure S3: Scatter plot of DHR values retrieved from CGLS, MODIS and MISR between 2012 and 2016. Figure S4: Scatter plot of BHR values retrieved from CGLS, MODIS and MISR between 2012 and 2016. Figure S5: Density plot of DHR and BHR upscaled from tower FoV to CGLS resolution over a 20\*20 pixel region.

**Author Contributions:** Conceptualization, J.-P.M. (Jan-Peter Muller); methodology, R.S. (Rui Song), S.K. (Said Kharbouche); software, R.S., S.K.; validation, R.S.; formal analysis, R.S.; investigation, R.S.; resources, J.-P.M.; data curation, J.-P.M.; writing—original draft preparation, R.S.; writing—review and editing, J.-P.M., S.K., and W.W. (William Woodgate); visualization, R.S.; supervision, J.-P.M.; project administration, J.-P.M.; funding acquisition, J.-P.M.

**Funding:** This research was funded by the European Commission Joint Research Centre grant number [FWC932059], part of the Global Component of the European Union's Copernicus Land Monitoring Service.

*Remote Sens.* **2019**, *11*, 644

**Acknowledgments:** This work used JASMIN, the UK's collaborative data analysis environment http://jasmin.ac. uk. This work has been undertaken using data from the Global Component of the European Union's Copernicus Land Monitoring Service, European Commission Joint Research Centre FWC [932059]. We would like to thank Nadine Gobron and Christian Lanconelli of JRC Ispra for fruitful discussions. We would like to thank NOAA for access to their datasets through SURFRAD (http://www.esrl.noaa.gov/gmd/grad/surfrad/) and the BSRN (http://bsrn.awi.de) [5] for access to their datasets. This work also used tower albedometer data acquired by the FLUXNET community and in particular by the following networks: AmeriFlux (U.S. Department of Energy, Biological and Environmental Research, Terrestrial Carbon Program (DE-FG02-04ER63917 and DE-FG02-04ER63911), CarboEuropeIP, CarboItaly, CarboMont, OzFlux, USCCC. This study has been undertaken using data from GBOV "Ground Based Observation for Validation" (https://land.copernicus.eu/global/gbov), part of the Global Component of the European Union's Copernicus Land Monitoring Service. GBOV product developments are managed by ACRI-ST from the research work of University College London, University of Leicester, University of Southampton, University of Valencia and Informus GmbH. We acknowledge the financial support to the tower albedometer data harmonization provided by CarboEuropeIP, FAO-GTOS-TCO, iLEAPS, Max Planck Institute for Biogeochemistry, National Science Foundation, University of Tuscia, Université Laval and Environment Canada and US Department of Energy and the database development and technical support from Berkeley Water Center, Lawrence Berkeley National Laboratory, Microsoft Research eScience, Oak Ridge National Laboratory, University of California-Berkeley, University of Virginia.

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


**4. Growth Rate Dynamics**

In this section, we make one modification to the delay best reply dynamical system, (1), and pursue the possibility of bounded dynamics when the system includes some nonlinearities. In particular, the growth rate adjustment is assumed in which the growth rate of output is controlled by the difference between the optimal output and the actual output,

$$\frac{\dot{x}\_i(t)}{x\_i(t)} = k \left[ \alpha - x\_i(t - \tau\_1) - \beta \sum\_{j \neq i}^{n} x\_j(t - \tau\_2) \right] \text{ for } i = 1, 2, \dots, n. \tag{11}$$

System (11) has the same stationary point as system (1). The homogeneous part of its linearized version is

$$\mathbf{x}\_{i}(t) = K \left[ -\mathbf{x}\_{i}(t - \tau\_{1}) - \beta \sum\_{j \neq i}^{n} \mathbf{x}\_{j}(t - \tau\_{2}) \right] \text{ for } i = 1, 2, \dots, n,\tag{12}$$

where

$$\mathbb{K} = kx^{\epsilon}.$$

Comparing (12) with (2) reveals that only the adjustment parameters are different. Thus, the formulas for the critical delays in (4), (5), (7) and (8) obtained in the best reply dynamic system can be applied to the growth rate dynamical system (12) if *k* is replaced with *K*.

The remaining part of this section is divided into two. The stability switching curves under the growth rate dynamics are constructed and numerical simulations are performed in the first subsection. The stability index is examined to provide theoretical backgrounds with the directions of stability switches for the numerical results in the second part.

## *4.1. Stability Switching Curves*

It is assumed henceforth that *K* replaces *k*. Then the pairs of (*τ*+ 1 (*ω*, -1), *τ*− 2 (*ω*, -2)) and (*τ*− 1 (*ω*, -1), *τ*+ 2 (*ω*, -2)) in (4) and (5) satisfy the following characteristic equation,

$$
\lambda + \mathcal{K}e^{-\lambda \tau\_1} - \mathcal{K}\beta e^{-\lambda \tau\_2} = 0 \tag{13}
$$

where the definitions of *θ*1 and *θ*2 should be changed to

$$\theta\_1(\omega) = \cos^{-1}\left[\frac{4\omega^2 + 3K^2}{8K\omega}\right],\tag{14}$$

$$\theta\_2(\omega) = \cos^{-1}\left[\frac{4\omega^2 - 3K^2}{4K\omega}\right] \tag{15}$$

and the interval *ω* is redefined by

$$I = \left[\frac{1}{2}K, \frac{3}{2}K\right].$$

We then have two sets of line segments in the first quadrant of the (*τ*1, *τ*2) plane,

$$L\_1^+(\ell\_1, \ell\_2) = \left\{ \left( \pi\_1^+(\omega, \ell\_1), \pi\_2^-(\omega, \ell\_2) \right) \mid \omega \in I, \ (\ell\_1, \ell\_2) \in \mathbf{Z} \right\} \tag{16}$$

and

*L*− 1 (-1, -2) = (*τ*− 1 (*ω*, -1), *τ*+ 2 (*ω*, -2)) *ω* ∈ *I*, (-1, -2) ∈ *Z* (17)

similar to the case of best reply dynamics. Lemma 1 characterizes the relations of the segments *L*+ 1 (-1, -2) and *L*− 1 (-1, -2) for the extreme values of *ω* in interval *I*.

**Lemma 1.** *L*+ 1 (-1, -2 +1) = *L*− 1 (-1, -2) *holds for the initial point of I, ω* = *K*/2, *and L*− 1 (-1, -2) = *L*+ 1 (-1, -2) *holds for the terminal point of I, ω* = 3*K*/2*.*

**Proof.** Substituting *ω* = *K*/2 into (14) and (15) gives

$$\theta\_1(\mathsf{K}/2) = \cos^{-1}(1) = 0 \text{ and } \theta\_2(\mathsf{K}/2) = \cos^{-1}(-1) = \pi\_1$$

implying that

$$
\pi\_1^{\pm}(\mathbb{K}/2, \ell\_1) = \frac{2}{\mathcal{K}} \left( \frac{3}{2}\pi + (2\ell\_1 - 1)\pi \right),
$$

and

$$
\pi\_2^+ (\mathsf{K}/2, \ell\_2 + 1) = \pi\_2^- (\mathsf{K}/2, \ell\_2) = \frac{2}{\mathsf{K}} \left( \frac{1}{2}\pi + 2\ell\_2\pi \right).
$$

Hence *L*+ 1 (-1, -2 + 1) = *L*− 1 (-1, -2) at the initial point of *I*. In the same way, for *ω* = 3*K*/2,

$$\theta\_1(3K/2) = \cos^{-1}(1) = 0 \text{ and } \theta\_2(3K/2) = \cos^{-1}(1) = 0$$

implying that

$$
\pi\_1^+(3K/2, \ell\_1) = \pi\_1^-(3K/2, \ell\_1) = \frac{2}{3K} \left(\frac{3}{2}\pi + (2\ell\_1 - 1)\pi\right).
$$

and

$$
\pi\_2^-\left(3K/2,\ell\_2\right) = \tau\_2^+\left(3K/2,\ell\_2\right) = \frac{2}{3K}\left(\frac{1}{2}\pi + \left(2\ell\_2 - 1\right)\pi\right).
$$

Hence *L*+ 1 (-1, -2) = *L*− 1 (-1, -2) at the terminal point of *I*. This completes the proof.

Pairs of (*τ*¯+ (*m*1), *τ*¯− (*m*1)) and (*τ*¯− (*m*2), *τ*¯+ (*m*2)) from (7) and (8) satisfy the characteristic equation,

$$
\lambda + K e^{-\lambda \tau\_1} + K \beta (n - 1) e^{-\lambda \tau\_2} = 0 \tag{18}
$$

where the definitions of ¯ *θ*1 and ¯ *θ*2 should be changed to

$$\theta\_1(\omega) = \cos^{-1}\left[\frac{4\omega^2 - K^2 \left(n - 3\right)\left(n + 1\right)}{8K\omega}\right] \tag{19}$$

and

$$\bar{\theta}\_2(\omega) = \cos^{-1}\left[\frac{4\omega^2 + K^2 \left(n - 3\right)\left(n + 1\right)}{4K(n - 1)\omega}\right] \tag{20}$$

and the interval for *ω* is defined, respectively, by

$$I\_2 = \left[\frac{1}{2}K, \frac{3}{2}K\right] \text{ if } n=2.$$

and

$$I\_{\mathbb{H}} = \left[\frac{n-3}{2}K, \frac{n+1}{2}K\right] \text{ if } n \ge 3.$$

We also have two line segments of (*τ*1, *τ*2),

$$L\_2^+\left(\mathcal{m}\_1, \mathcal{m}\_2\right) = \left\{ \left(\overline{\mathfrak{r}}\_1^+\left(\omega, \mathfrak{m}\_1\right), \overline{\mathfrak{r}}\_2^-\left(\omega, \mathfrak{m}\_2\right)\right) \mid \omega \in I\_2 \text{ or } I\_{\mathfrak{n}\_1}\left(\mathfrak{m}\_1, \mathfrak{m}\_2\right) \in \mathbf{Z} \right\} \tag{21}$$

and

$$L\_2^-\left(m\_1, m\_2\right) = \left\{ \left(\mathfrak{r}\_1^-\left(\omega, m\_1\right), \mathfrak{r}\_2^+\left(\omega, m\_2\right)\right) \mid \omega \in I \text{ or } I\_{n\prime}\left(m\_1, m\_2\right) \in \mathbf{Z} \right\} \tag{22}$$

similarly to the case of best reply dynamics. Similarly to Lemma 1, we have the followings:

**Lemma 2.** *In the case of n* = 2, *L*+ (*m*1, *m*2 + 1) = *L*− (*m*1, *m*2) *holds for the initial point of I*2*, ω* = *K*/2, *and L*− (*m*1, *m*2) = *L*+ (*m*1, *m*2) *holds for the terminal point of I*2*, ω* = 3*K*/2*.*

Notice that for *n* = 3,

$$\lim\_{\omega \to 0} \tau\_1^{\pm}(\omega, m\_1) = \infty \text{ and } \lim\_{\omega \to 0} \tau\_2^{\pm}(\omega, m\_2) = \infty.$$

The equality of the segments does not hold at the initial point of *I*3 but only at the terminal point which can be proved similarly to Lemma 1.

**Lemma 3.** *In the case of n* = 3, *L*+ 2 (*m*1, *m*2) = *L*− 2 (*m*1, *m*2) *holds for the terminal point of I*3*, ω* = 2*K*.

If *n* ≥ 4, then the following result holds.

**Lemma 4.** *In the case of n* ≥ 4, *L*+ (*m*1, *m*2) = *L*− (*m*1 + 1, *m*2) *holds for the initial point of In, ω* = (*n* − 1)*K*/2 *and L*+ (*m*1, *m*2) = *L*− (*m*1, *m*2) *holds for the terminal point ω* = (*n* + 1)*K*/2*.*

In the following, we will construct stability switching curves. To this end, we specify the parameters' values as *α* = 10 and *k* = 0.1. In Figure 3, the dotted red loci are described by *L*− 1 (-1, -2) with -1 = 0 and -2 = 0, 1 and the dotted blue locus by *L*+ 1 (0, 1). The black point *a* is the initial point of *L*+ 1 (0, 1) and *L*− 1 (0, 0) and its coordinates are

$$a' = \left(\frac{3}{2}\pi, \frac{3}{2}\pi\right).$$

at which *L*+ 1 (0, 1) = *L*− 1 (0, 0) holds by Lemma 1. The black point *b* is the terminal point of *L*+ 1 (0, 1) and *L*− 1 (0, 1) and its coordinates are

$$b' = \left(\frac{1}{2}\pi, \frac{3}{2}\pi\right).$$

at which *L*+ 1 (0, 1) = *L*− 1 (0, 1) holds by Lemma 1. The blue and red solid curves are described by *L*+ 2 (0, 0) and *L*− 2 (0, 0). They are connected at point *a*,

$$a = \left(\frac{1}{2}\pi, \frac{1}{2}\pi\right).$$

at which *L*+ 2 (0, 0) = *L*− 2 (0, 0) by Lemma 2.

The dotted and solid curves are smoothly connected as is seen in Figure 3. As a result, the (*τ*1, *τ*2) region is divided into two subregions by the stability switching curve connecting the left-most parts among the segments of *L*± 2 (0, 0), *L*± 1 (0, 1), and *L*− 1 (0, 1). As the Cournot equilibrium is stable when there are no delays, it is stable in the region including the origin and left to the connecting curve.

**Figure 3.** Stability switching curve (SSC) with *n* = 2.

We want to investigate the influence of *τ*1 and *τ*2. Two simulations in the case of *n* = 2 are performed with initial functions,

$$
\varphi\_1(t) = \mathfrak{x}\_1^\varepsilon - 2 \text{ and } \varphi\_2(t) = \mathfrak{x}\_2^\varepsilon + 1 \text{ for } t \le 0.
$$

The first simulation result along the diagonal is presented in Figure 4A. The delays increase from *τi* = 0 to *τi* = 3.4 with an increment of 0.003 along the diagonal. The Cournot equilibrium is asymptotically stable for smaller delays and becomes unstable through a Hopf bifurcation at

$$
\tau\_i^a = \tau\_0^\*(2) = \frac{1}{2}\pi \text{ for } i = 1, 2, 3
$$

producing a limit cycle that further bifurcates to a multi-periodic cycle for larger delays. The second result with the different two delays is given in Figure 4B. The value of *τ*1 increases from *τ*1 = *τA* 1 ( 1.423) to *τ*1 = *τB* 1 (= 3.4) along the dotted horizontal line at *τ*2 = 3. More precisely, the bifurcation diagrams with two delays are constructed in the following procedure with *Mathematica*, version 12.1. The value of *τ*2 is fixed at 3, and the value of *τ*1 is increased from *τ*min = *τA* 1 to *τ*max = *τB* 1 with an increment (*τ*max − *τ*min)/1000. For each value of *τ*1, dynamic system (11) runs for 0 ≤ *t* ≤ *T*(= 2000), and the data for *t* ≤ *T* − 100 are discarded to get rid of the initial disturbance. The local maxima and minima out of the remaining data are plotted against this *τ*1 value. Then the value of *τ*1 is increased and then the same procedure is repeated until *τ*1 arrives at *τ*max. The following bifurcation diagrams are obtained in the same way. The resulting bifurcation diagram shows that the dynamic system experience similar dynamics. The stability of the equilibrium point is confirmed for the zero delay and holds for *τ*1 < *τA* 1 and *τ*2 = 3. In both diagrams (and the following diagrams), notation *x*˜*e* = log [*xe* ] is used.

**Figure 4.** Bifurcation diagrams with *n* = 2.

We now increase the number of firms to 3. Figure 5A shows the stability switching curves. The line segments of *L*+ 2 (0, 0) (i.e., the solid blue curve) and *L*− 2 (0, 0) (i.e., the solid red curve) take the *L*-shaped profile and rotate counter-clockwise at point *a* to the extent that the solid red curve is located furthermost to the left. By Lemma 4, both line segments head to point *a*, the terminal point as *ω* increases to 2*K*. We simulate the model (11) along the diagonal (i.e., *τ*1 = *τ*2) and the dotted horizontal line at *τ*2 = 3 (i.e., *τ*1 = *τ*2) in Figure 5A. As we find qualitatively no big differences between these simulation results as in Figure 4A,B, we depict only the bifurcation diagram with different delays in Figure 5B. It is seen that alá "period-doubling bifurcation" occurs in which the Cournot equilibrium is

asymptotically stable for *τ*1 < *τA* 1 (1.136), loses stability at *τ*1 = *τA* 1 and bifurcates to a limit cycle from which new limit cycles emerge having a doubled period of the cycle as *τ*1 increases from *τA* 1 . We also see that further increasing *τ*1 gives rise to complicated dynamics that suddenly shrinks to a limit cycle with multiple local maxima and minima at some critical point.

**Figure 5.** Dynamic properties of Equation (11) with *n* = 3.

In the case of *n* = 4, as is seen in Figure 6A, the solid red and blue segments rotate counter-clockwise further at point *a*, leading to that the red segment crosses the vertical axis. In Figure 6B, we see that the bifurcation diagram gets more complicated and various dynamics can emerge.

**Figure 6.** Dynamic properties of Equation (11) with *n* = 4.

Lastly, we simulate system (11) with *n* = 9. The shape of the stability switching curve is different from those with smaller *n*. In Figure 7A, the positive-sloping dotted line is the diagonal, the dotted-red line is *L*− 2 (0, 0) as before and the black dots are the starting or ending points of the segments. A remarkable difference is that the solid red-blue segments consist of the wave-shaped curve. Accordingly, the bifurcation diagram is obtained along the horizontal dotted line at *τ*2 = 2 and exhibits a different route to chaos. The stability of the Cournot equilibrium is lost at *τ*1 = *τA* 1

(0.646), regained at *τ*1 = *τB* 1 (5.441), and then lost again at *τ*1 = *τC* 1 (7.306). Unstable oscillatory trajectories get complicated for *τ*1 > *τD* 1 (7.697). It is known that time delays destabilize dynamic systems. This simulation, however, indicates that time delays can also stabilize the systems.

**Figure 7.** Dynamic properties of Equation (11) with *n* = 9.

Dynamic system (11) examines the birth of complicated dynamics through a period-doubling bifurcation and the occurrence of stability loss and gain. Needless to say, time delays play prominent roles. In addition, taking account of the fact that only the firm's number is different in those numerical studies, the larger number could influence the system's dynamics by increasing the degree of interactions among the firms.

## *4.2. Stability Index*

We compute the stability index to provide a theoretical background for finding directions of stability switches. First, we denote the second and third vectors of (3) by *Q*1 and *Q*2,

$$Q\_1 = a\_1(i\omega)e^{-i\omega\tau\_1} = -i\frac{K}{\omega} \left(\cos\omega\tau\_1 - i\sin\omega\tau\_1\right).$$

and

$$Q\_2 = a\_2(i\omega)e^{-i\omega\tau\_2} = i\frac{\mathbb{K}\beta}{\omega} \left(\cos\omega\tau\_2 - i\sin\omega\tau\_2\right).$$

Having *Q*1 and *Q*2, we further denote the real and imaginary parts by the followings:

$$R\_1 = \text{Re}Q\_1 = -\frac{K}{\omega}\sin\omega\,\tau\_1 \text{ and } I\_1 = \text{Im}Q\_1 = -\frac{K}{\omega}\cos\omega\,\tau\_1$$

$$\text{and}$$

$$R\_2 = \text{Re}Q\_2 = \frac{\text{K}\beta}{\omega} \sin \omega \tau\_2 \text{ and } I\_2 = \text{Im}Q\_2 = \frac{\text{K}\beta}{\omega} \cos \omega \tau\_2.$$

Finally, the stability index is defined as follows:

$$\begin{array}{rcl} S &=& R\_2 I\_1 - R\_1 I\_2 \\ &=& \frac{K^2 \beta}{\omega^2} \left( \sin \omega \tau\_1 \cos \omega \tau\_2 - \cos \omega \tau\_1 \sin \omega \tau\_2 \right) . \end{array}$$

hence

$$S = \frac{K^2 \beta}{\omega^2} \sin\left[\omega \left(\tau\_1 - \tau\_2\right)\right]. \tag{23}$$

In the same way, we denote the second and third vectors of (6) by *Q*¯ 1 and *Q*¯ 2,

$$Q\_1 = b\_1(i\omega)e^{-i\omega\tau\_1} = -i\frac{K}{\omega} \left(\cos\omega\tau\_1 - i\sin\omega\tau\_1\right).$$

and

$$\bar{Q}\_2 = b\_2(i\omega)e^{-i\omega\tau\_2} = -i\frac{K\beta(n-1)}{\omega} \left(\cos\omega\tau\_2 - i\sin\omega\tau\_2\right).$$

The real and imaginary parts are the followings:

$$\mathcal{R}\_1 = \text{Re}\mathcal{Q}\_1 = -\frac{K}{\omega}\sin\omega\,\text{tr}\_1 \text{ and } I\_1 = \text{Im}\mathcal{Q}\_1 = -\frac{K}{\omega}\cos\omega\,\text{tr}\_1$$

and

$$\bar{R}\_2 = \text{Re}\bar{Q}\_2 = -\frac{K\beta(n-1)}{\omega}\sin\omega\tau\_2 \text{ and } \bar{I}\_2 = \text{Im}\bar{Q}\_2 = -\frac{K\beta(n-1)}{\omega}\cos\omega\tau\_2.$$

moreover, the stability index is as follows:

$$\begin{aligned} \vec{S} &= -\mathcal{R}\_2 I\_1 - \mathcal{R}\_1 I\_2 \\ &= -\frac{K^2 \beta (n-1)}{\omega^2} \left( \sin \omega \tau\_1 \cos \omega \tau\_2 - \cos \omega \tau\_1 \sin \omega \tau\_2 \right) . \end{aligned}$$

Hence

$$S = -\frac{K^2 \beta (n - 1)}{\omega^2} \sin \left[ \omega \left( \tau\_1 - \tau\_2 \right) \right]. \tag{24}$$

We call the direction of the curve that corresponds to increasing *ω* the *positive direction*. We also call the region on the left-hand side *the region on the left* when we head in the positive direction of the curve. *Region on the right* is defined similarly. Concerning the stability changes, we have the following result from Matsumoto and Szidarovszky (2018) that is based on Gu et al. (2005):

**Theorem 2.** *Let* (*τ*1, *τ*2) *be a point on the stability switching curves, when iω is a simple pure complex eigenvalue. Assume we look toward increasing values of ω on the curve, and a point* (*τ*1, *τ*2) *moves from the region on the right to the region on the left. A pair of eigenvalues crosses the imaginary axis to the right if S* > 0 *or S*¯ > 0*. If S* < 0 *or S*¯ < 0*, then crossing is in the opposite direction.*

The condition of the theorem is satisfied if all *iω* egenvalues are single. It can be proved that the multiple eigenvalues, if any, are isolated from each other, so do the corresponding points on the stability switching curve. Hence at these points, the directions of stability switching are the same as those in the points of their neighborhoods.

We now compute the stability index on the solid red segment of the stability switching curve in Figure 3. The red segment is a locus of the following points,

$$L\_2^-(0,0) = \left\{ \left. \left( \mathfrak{t}\_1^-(\omega,0), \mathfrak{t}\_2^+(\omega,0) \right) \right| \: \omega \in \left[ \frac{1}{2} K, \frac{3}{2} K \right] \right\}.$$

From (7) and (8), we have

$$\begin{aligned} \omega \left( \bar{\tau}\_1^- (\omega, 0) - \bar{\tau}\_2^+ (\omega, 0) \right) &= \begin{bmatrix} 3\\ 2\pi - \pi - \bar{\theta}\_1 (\omega) \end{bmatrix} - \begin{bmatrix} 3\pi\\ 2 \end{bmatrix} - \pi + \bar{\theta}\_2 (\omega) \begin{bmatrix} 1\\ 2 \end{bmatrix} \\ &= \begin{bmatrix} \theta\_1 (\omega) + \bar{\theta}\_2 (\omega) \end{bmatrix} \end{aligned}$$

implying

$$\sin\left[\omega\left(\pi\_1-\pi\_2\right)\right] = -\sin\left[\bar{\theta}\_1(\omega) + \bar{\theta}\_2(\omega)\right] < 0$$

when *θ*1 + *θ*2 < *π*. If *θ*1 + *θ*2 = *π*, then the triangle reduces to a line such that

$$|a\_1(i\omega)| - |a\_2(i\omega)| = \pm 1.$$

That is, in Equation (13),

$$\frac{K}{\omega} - \frac{K\beta}{\omega} = \frac{K}{2\omega} = 1$$

showing that *ω* = *π*/2 being the left endpoint of interval *I*, given for *ω*, which gives the common starting point of two line segments. In the case of Equation (18),

$$\frac{K}{\omega} - \frac{K\beta(n-1)}{\omega} = \frac{K}{\omega} \left(\frac{3-n}{2}\right).$$

If *n* = 2, this equals +1 if *ω* = *K*/2, which is the initial point of *I*. If *n* = 3, then this expression is always zero, so cannot be +1 or −1. If *n* > 3, then this expression can be only −1, when *ω* = *K*(*n* − 3)/2, which is the left endpoint of interval *In* which gives again the common starting point of two line segments. In these points, the direction of stability switching is the same as that in the two connecting segments. So in the rest of the discussion, we will assume that *θ*1 + *θ*2 < *π*. Hence the stability index *S*¯ is positive on the solid red segments of the stability switching curve. In Figure 3, the arrows on the solid red segment indicate the positive direction and the red *R* and *L* mean the right and left regions along the red segment. As (*τ*1, *τ*2) moves from the *R*-region to the *L*-region and *S*¯ > 0, Theorem 2 implies that a solution pair of (18) crosses the imaginary axis to the right. That is, stability is lost. As seen in Figure 4B, the stability is lost at point *A* with *τ*1 = *τA* 1 when *τ*1 increases along the horizontal dotted line at *τ*0 2 = 3.

Similarly, we can compute the stability index on the solid blue segment,

$$L\_2^+(0,0) = \left\{ \left( \vec{\tau}\_1^+(\omega,0), \vec{\tau}\_2^-(\omega,0) \right) \mid \omega \in \left[ \frac{1}{2}K, \frac{3}{2}K \right] \right\}.$$

From (7) and (8) with *K*,

$$\begin{aligned} \omega \left( \mathfrak{t}\_1^+ (\omega, 0) - \mathfrak{t}\_2^- (\omega, 0) \right) &= \begin{bmatrix} \frac{3}{2} \pi - \pi + \bar{\theta}\_1 (\omega) \\ \end{bmatrix} - \begin{bmatrix} 3\pi \\ 2 \end{bmatrix} - \pi - \bar{\theta}\_2 (\omega) \begin{bmatrix} \\ \end{bmatrix} \\ &= \begin{array}{c} \bar{\theta}\_1 (\omega) + \bar{\theta}\_2 (\omega) . \end{bmatrix} \end{aligned}$$

Then

$$\sin\left[\omega\left(\bar{\tau}\_1^+(\omega,0)-\bar{\tau}\_2^-(\omega,0)\right)\right] = \sin\left[\bar{\theta}\_1(\omega)+\bar{\theta}\_2(\omega)\right] > 0.$$

The stability index *S*¯ is negative,

$$\bar{S} = -\frac{K^2 \beta (n - 1)}{\omega^2} \sin \left[ \bar{\theta}\_1(\omega) + \bar{\theta}\_2(\omega) \right] < 0.1$$

The blue *L* and *R* denote the right-region and the left-region with respect to the solid blue segment. Hence the stability is lost when a pair of (*τ*1, *τ*2) crosses the blue segment from the *L*-region to the *R*-region.

Consider the stability switching on the dotted red segment located in the upper-left corner of Figure 3. The segment is described by

$$L\_1^-(0,1) = \left\{ \left( \pi\_1^-(\omega,0), \pi\_2^+(\omega,1) \right) \mid \omega \in \left[ \frac{1}{2}K, \frac{3}{2}K \right] \right\}.$$

Then

$$
\omega \left( \mathfrak{r}\_1^- (\omega, 0) - \mathfrak{r}\_2^+ (\omega, 1) \right) = -\pi - \left( \theta\_1(\omega) + \theta\_2(\omega) \right).
$$

The stability index is positive

$$S = \frac{K^2 \beta}{\omega} \sin \left[\theta\_1(\omega) + \theta\_2(\omega)\right] > 0$$

showing that crossing these segments from *R* to *L*, stability is lost.

In the lower part of Figure 3, there is a small segment of *L*− 1 (0, 0) where

$$L\_1^-(0,0) = \left\{ \tau\_1^-(\omega,0), \ \tau\_2^+(\omega,0) \mid \omega \in \left[\frac{1}{2}K, \frac{3}{2}K\right] \right\},$$

so

$$\begin{aligned} \omega \left( \tau\_1^- (\omega, 0) - \tau\_2^+ (\omega, 0) \right) &= \left[ \frac{3}{2} \pi - \pi - \theta\_1 (\omega) \right] - \left[ \frac{1}{2} \pi + \pi - \theta\_2 (\omega) \right] \\ &= \left[ \pi - \left[ \theta\_1 (\omega) + \theta\_2 (\omega) \right] \right] .\end{aligned}$$

Then

$$\sin\left[\omega\left(\pi\_1^-\left(\omega,0\right)-\pi\_2^+\left(\omega,0\right)\right)\right] = \sin\left[\bar{\theta}\_1(\omega)+\bar{\theta}\_2(\omega)\right]>0$$

meaning that crossing this segment from the stable region, at least one eigenvalue changes the sign of its real part from negative to positive, implying stability loss.

#### **5. Concluding Remarks**

In this paper, *n*-firm dynamic oligopolies were examined without product differentiation and with linear price and cost functions. Continuous time scales were assumed reconsidering the classical dynamic model of McManus and Quandt (1961) with the best response dynamics. Without delays, the equilibrium is always asymptotically stable without delays regardless of the values of the positive adjustment speeds. We examined how this stability is lost when the firms face implementation and information delays. For the sake of mathematical simplicity, it was assumed that the firms have the same marginal costs and identical delays in both types. If these delays are equal, then a single-delay model is obtained. If the delay is sufficiently small, then the equilibrium is oscillatory stable, at the threshold, the trajectories show cyclic behavior and for larger delays, the cycles become expanding. If the delays are different, then in the resulting two-delay case the stability switching curves were first constructed and then the directions of the stability switches were determined. Growth rate dynamics brought nonlinearities into the model, but their linearized version is identical with best response dynamics, so shows similar local dynamics. Numerical results and simulation studies verify and illustrate the theoretical findings.

This research can be continued in two different ways. One is the consideration of different model modifications such as product differentiation, multi-product models, oligopsonies, labor-managed, and rent seeking oligopolies, including market saturation to mention only a few. The other research direction could be to examine nonlinear models, the local dynamics are similar to that of linear models, however with very different global dynamic behavior.

**Author Contributions:** Conceptualization, methodology, A.M. and F.S.; software, A.M.; validation, A.M. and F.S.; formal analysis, A.M.; writing—original draft preparation, A.M.; writing—review and editing, F.S.; visualization, A.M.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received the financial support from the Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research (C) 20K01566).

**Acknowledgments:** The authors would like to thank the anonymous reviewers for their careful reading and valuable comments.

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


**3.4. Dandelion (***Taraxacum officinale***)**

Dandelion is a perennial herb of multiple global varieties that has traditionally been used for liver, spleen, kidney, and gastrointestinal disorders, although there have been no clinical trials investigatingtheeffectsofdandelioninweightloss [27,178].Itiscommonlyaddedtoweightloss supplements,mainlyforitsdiureticproperties,althoughtheherbdoespossesssomemildlaxative properties [179-181].Therearenoknowndrugs interactionsbetween*Taraxacum*andotherherbs or drugs, although one study in rats suggests a probable interaction with quinolone antibiotics due to the high mineral content of *Taraxacum* [182]. In the study, ciprofloxacin (20 mg/kg) Cmax significantly decreased while Vd and t½ significantly increased when administered with crude dandelion extract (2 g/kg) compared to control. There is one case report of hypoglycemia in a 58 year old diabetic patient following a 2-week period of dandelion consumption in salads [183]. Thepatientdeniedchanges incalorie consumption, exercise,orinsulindosing.Diabeticpatients taking hypoglycemic agents while consuming dandelion should be monitored.