Literature DB >> 34055583

Stability analysis and simulation of the novel Corornavirus mathematical model via the Caputo fractional-order derivative: A case study of Algeria.

Yacine El Hadj Moussa1, Ahmed Boudaoui2, Saif Ullah3, Fatma Bozkurt4, Thabet Abdeljawad5,6,7, Manar A Alqudah8.   

Abstract

The novel coronavirus infectious disease (or COVID-19) almost spread widely around the world and causes a huge panic in the human population. To explore the complex dynamics of this novel infection, several mathematical epidemic models have been adopted and simulated using the statistical data of COVID-19 in various regions. In this paper, we present a new nonlinear fractional order model in the Caputo sense to analyze and simulate the dynamics of this viral disease with a case study of Algeria. Initially, after the model formulation, we utilize the well-known least square approach to estimate the model parameters from the reported COVID-19 cases in Algeria for a selected period of time. We perform the existence and uniqueness of the model solution which are proved via the Picard-Lindelöf method. We further compute the basic reproduction numbers and equilibrium points, then we explore the local and global stability of both the disease-free equilibrium point and the endemic equilibrium point. Finally, numerical results and graphical simulation are given to demonstrate the impact of various model parameters and fractional order on the disease dynamics and control.
© 2021 The Authors.

Entities:  

Keywords:  COVID-19; Lyapunov function; Real data; Simulations; Stability analysis

Year:  2021        PMID: 34055583      PMCID: PMC8141347          DOI: 10.1016/j.rinp.2021.104324

Source DB:  PubMed          Journal:  Results Phys        ISSN: 2211-3797            Impact factor:   4.476


Introduction

Coronavirus is an infectious disease that affects the respiratory tracts and the strain is a virus known as SARS-Cov-2. The first case of this novel disease was reported in December 2019, in Wuhan, a city in China [1]. Factors such as hand contacts, sneezes, talk and breaths of an infectious person are main causes to spread the virus [2]. The novel COVID-19 is declared as a pandemic on 11 March 2020 by WHO and it affects almost all countries all over the world. It is considered to be a highly contagious disease affecting about 61 million infection cases and about 1.5 million death in a recent estimations. Although, the specific symptoms are not identified due to variability in its nature, however, the common symptoms appeared so far are fever and cough in addition to fatigue, difficult breathing etc. The COVID-19 is usually dangerous for aged people and for a person who has some Chronic diseases like tuberculosis or other respiratory problems. There is an incubation period of this infectious disease which is the moment of a person becomes infected for the first time and the appearance of the first symptoms. This period varies between three to six days [3]. In some cases, a person infected with the virus but do not develop symptoms at any moment in time, are classified as asymptomatic and they can spread the virus [4], [5]. Currently, no proper and effective treatment is available for a COVID-19 infected persons except for some drugs as Remdesivir which are approved by some countries like Australia and the European Union [6]. There is no effective and authorized vaccine for this novel infection although few countries have claimed it. The best prevention strategies used in almost countries to decrease and delay the epidemic pick (flattening the curve) are the frequent tests to determine the infected persons, isolation and lockdown, social distancing, use of strict SOPs, etc until effective treatments and vaccine become available. The preventive measures are only the way to reduces the chance of infections and slow the spread of the virus. The researcher around the globe are focusing to plain a strategy to overcome the COVID-19 pandemic. For this purpose, they have been adopted different approaches in order to explore the complex transmission dynamics of this infection. Mathematical models are one of considerable tools in this regard. Many epidemic models were introduced to analyze the dynamic of the COVID-19 and discuss the possible control strategies for the disease elimination in different regions of the world. For instance, a model of COVID-19 with Lockdown is suggested in [7], and the impact of undetected cases via a mathematical model is explored in [8]. The impact of some preventive measure on the curtailing of COVID-19 in Pakistan via a new mathematical model is presented in [9]. A transmission mathematical model considering the environmental spread of the virus with a case study of Saudi Arabia is studied in [10]. Most of the models of COVID-19 are formulated in terms of the integer order derivatives which have some limitations to describe the realistic aspects of a phenomena under consideration. In order to deal with those limitations, fractional order derivatives provide a useful tool for the disease dynamic and give some additional results that need to understand the models. Fractional order mathematical models possess memory propriety and provide a better scenario to describe an epidemic model. Many models on the dynamics different diseases in term of fractional order derivatives were proposed see for instance [11], [12], [13], [14], [15], [16], [17] and the literature referenced therein. Recently, many mathematical epidemic models with fractional order have been developed to investigated the dynamical patterns of the novel COVID-19 pandemic. For instance, in [18] a model of COVID-19 is formulated using Caputo fractional operators in which the authors predicted the situation of COVID-19 in the mainland of China. A fractional mathematical models based on the COVID-19 dynamics was proposed in [19] where the authors used Atangana-Baleanu operator to construct the model. In [20] a fractional order epidemic model is developed to describe the interaction among the bats and unknown hosts. In [21] the authors investigate the pandemic dynamics in India using SEIR compartmental model based on a non-singular derivative. A novel iterative scheme known as q- Homotopy analysis transform method is used to investigate the solution of proposed model. The dynamic and numerical approximation for the arbitrary-order coronavirus disease system is developed in [22] by considering three various arbitrary-order derivative operators. A study in [23] extend a classical model to fractional model and it is concluded that the model with non integer orders gives more interesting and deeper characteristics of the disease dynamics. In [24] the authors presented a robust numerical assessment of COVID-19 infection mathematical model by hermite wavelets scheme and they give a comparative analysis between solution obtained by Hermite wavelets and the solution by the ABM predictor corrector scheme. In [25] the authors explored the dynamical behavior of COVID-19 with some effective controlling measure to mitigate the infection. The purpose of the present article is to investigate the dynamics of the novel COVID-19 pandemic in Algeria using a compartmental model via the Caputo fractional order derivatives. To construct the model the infected population is divided into reported and unreported classes. The evolution of contagious person according the reported category which are detected by doing the test of COVID 19, and the unreported cases which concern people who have not performed the test. Furthermore, we make the hypothesis that reported cases migrate to critical cases before recovery or die. This approach is justified according the limitation to precede to testing all persons in some countries. The parameters of the proposed mathematical model for COVID-19 are estimated from the real data that was reported in Algeria. The rest of this paper is organized: The preliminary results are given in Section 2 The Caputo fractional order model and estimation procedure are given in Section 3. In Section 4, we deal with the basic analysis of the model such as the existence and uniqueness of the model solution in addition to the basic reproduction number which give information about the contagiousness or transmissibility of infectious diseases stability of the model equilibria etc. Then we deal in Section 5 with the numerical simulations and discussion of the model. The last section of the study presents the concluding remarks.

Preliminary results

We give basic definition and features of fractional calculus, for more details we can see [26]. The well-known Riemann Louiville integral with non integer order say of a function is presented bywhere , is the order of integration, and , is the Gamma function. The classical Caputo derivative with non integer order say , of a function is presented bywhere , with is the integer part of Caputo derivative and the Riemann Liouville integral satisfy the following properties , where , , If is such that , then .

Mathematical model description

In this section, we introduce a fractional order COVID-19 transmission model to investigate the dynamical behavior of outbreak of COVID-19 pandemic with respect to different values of parameters. We extend the SAIRD transmission model by introducing three infection sub-classes of I namely the undetected infected people, the detected infected people and critical infected people. Thus the total population shown by at time t is classified into variables and which denote various population sub-classes representing the number of susceptible, asymptomatic, undetected infected, detected infected, critical infected,recovered people from the virus and finally the number of dead people by the disease respectively. The well-known Caputo derivative with order is utilized to formulate the propped epidemic model that give the dynamic of the virus in different compartments summarized by a diagram presented in Fig. 1 . Thus the transmission model is organized by the following fractional differential equations
Fig. 1

scheme of different stages of transmission of a novel corona virus in different compartments.

scheme of different stages of transmission of a novel corona virus in different compartments. In addition, the following initial conditions are taken in consideration: In the model, a net inflow of susceptible people into the region is shown by at unit time. However, the susceptible people decreases after acquiring infection due to interaction with undetected infected people, detected infected people or critical infected people. A newly infected individual from compartment S becomes asymptomatic where the parameters , and are the transmission rates relevant undetected, detected and critical infected people respectively. The asymptomatic people progress through the undetected infectious and detected infectious compartments with an average , where a faction of asymptomatic people move to undetected infectious compartment, while a fraction of asymptomatic people move to detected infectious compartment. The reason for considering the classes and is due to the fact that we are not able to do massive COVID-19 test to the population. The reported infected individuals are based on the evolution of doing the PCR test and the unreported cases concerns people who have not performed the PCR test. Then detected infectious compartment will becomes in critical situation and move to infected critical compartment at a rate . The parameters and are the disease caused death rate for compartment and respectively. Notice that the present of the virus caused death in compartment is because in Algeria, all the dead people are tested COVID-19. Thus, even the undetected infected can be known that he died of the virus, also the present of the virus caused death in compartment is because we assume here that infected person in the compartment move to before they’re dead by COVID-19. Finally is the natural caused death in all compartment. Finally some of the undetected infectious people, detected infectious people and critical infectious people progress to recovery compartment at a rate and respectively, while some of them move to death virus compartment D at the rate respectively. Again notice that the progression of the undetected infectious people to recovery compartment is due to some person who develop self immunity for the virus. The last equation of the model (3) denotes the total death cases due to COVID-19.

Parameter estimations procedure

This part of the paper present the estimation procedure of the model parameters. For this purpose, a standard least squares approach is utilized. This approach is used to find the best fit for a set of data points by minimizing the sum of the offsets or residuals of points from the plotted curve and has been applied in many recent literature [9], [12]. The birth rate and the natural mortality rate denoted by were estimated from the literature [27]. The model fitting curve to the real data points is shown in Fig. 2 while the resulting estimated as well as the fitted parameters are listed in Table 1 . The numerical value of the most important threshold parameter known as the basic reproduction number is approximated as 1.20.
Fig. 2

Model predicted infected curve (solid blue line) to the confirmed infected cases in Algeria with the parameter .

Table 1

Model parameters with biological meaning and respective fitted values.

Parameter symbolsBiological meaningValuereference
ΔRecruitment rate1534Estimated[27]
μNatural mortality rate1/(77.5*365)[27]
ν1Transmission rate of undetected people0.5022Fitted
ν2Transmission rate of detected people0.4666Fitted
ν3Transmission rate of critical people0.4080Fitted
d1Death rate in Iu class due to infection0.0575Fitted
d2Death rate in Ir class due to infection0.0271Fitted
d3Death rate in Ic class due to infection0.0141Fitted
σIncubation rate0.6151Fitted
γIuRecovery rate of undetected infectious people0.3333Fitted
γIrRecovery rate of detected infectious people0.6026Fitted
γIcRecovery rate of critical infectious people0.3103Fitted
δIrCritical rate of detected infectious people0.4972Fitted
Model predicted infected curve (solid blue line) to the confirmed infected cases in Algeria with the parameter . Model parameters with biological meaning and respective fitted values.

Mathematical analysis of the model

Adding up the equations given in (3), we have as . The above inequality leads to Before showing the existence and uniqueness of solution to the system Eqs. (3), (4), note that the equation involving D does not appear in the other equations, and since the population is closed we can determine D by, So, we restrict our self throughout the paper to the model with the three subclass of I namely , and . Hence we consider the following fractional differential equationsand the following initial conditions: Let , and . We have the following result about the solution and the feasible region for the model described in (3). There is a unique solution for the initial value problem given by (3), (4) , and the solution hold for all in Now, to show the existence and uniqueness we use fixed point theory and Picard-Lindelöf method. To proceed, we may rewrite the system described in (8) in the following classical formwhere, the vector and the function is defined as followsand the initial condition is To do so we proceed in the following manner. Using initial conditions (12) and fractional integral operator 4, we transform the system (8) into the following integral equations Using (11) and the definition of in (1), we obtained the state variable in terms of , where . The Picard iterations are given by the following equations Corresponding to the form (10), and with the initial condition (12) we have the following integral equation The function defined in (13) satisfies the Lipschitz condition given by where and the norm corresponds to the space . We proof lemma only for the other can be obtained in the same manner. Then, Also we have  □ Assuming we have (14) , then there exist a unique solution to the system (8), (9) if The solution to the system (8), (9) iswhere, T is the Picard operator defined by Further, we have Since, then, the operator T is a contraction, hence the system (8), (9) has a unique solution. □ For the non negative solution we need the following lemma.

Generalized mean value theorem

Let and for , then we have with , for all . From Lemma 4.4 we have If then the function f is non decreasing for all . If then the function f is non increasing for all .

Proof of Theorem 4.1

For the existence and uniqueness we use Lemma 4.2 and Lemma 4.3. Now to proof that solutions of the system (8), (9) with non negative initial data will remain non negative for all , we use Lemma 4.4. Since Then we conclude that the solution of system (8), (9) belongs to . □

Equilibrium points

To evaluate the disease free equilibrium (DFE) let. . System (8) become For which we get the DFE as follows:

Basic reproduction number and interpretation

In epidemiology the expected number is the basic reproduction number denoted by , it is useful to determine if the infection can spread in a population, witch lead to the strategies for eradicate the infection. One approach used to evaluate is the next generation matrix introduced in [28]. Let , the system (8) can be written aswhere, The respective Jacobian of above matrices at the DFE are evaluated as follows: Hence, from the result in [28], the basic reproduction number after some calculations is obtain as:

Interpretation of

The basic reproductive number is the summation of three terms i.e., . The first term gives the number of new infection cases moved from asymptotic compartment A to the class. Clearly, it composed of the product of infection transmission rate in compartment denoted by , the fraction of asymptomatic cases that completed the incubation period and migrate to undetected person with the average period spend in A compartment and the average period spend in the compartment . The second term represents the number of new COVID-19 infections moved from asymptotic compartment A to the compartment, and similarly the third term express the number of new COVID-19 infections moved from asymptotic compartment A to the class compartment.

Stability analysis

The following results provide the local and global stability results of the system (8) around the DFE . The DFE of the system (8) is locally asymptotically stable when . The associated Jacobian matrix of (8) evaluated at is given bywhere Computations give the following characteristic polynomial In the characteristic polynomial of we have the eigenvalues Which show negative real parts if . This complete the proof. □ With the notation in (27) system (8) becomesand in addition the number become We establish result of global stability for fractional differential system (8). For the proof of the Theorem 4.7 we need this lemma: [Lemma 2.4 [29]] Let denotes a continuous and derivable function then, for any time instant . The dDFE for the system (8) is globally asymptotically stable, if . Let and consider the following Lyapunov functionwhere and are defined in (27). Using linearity of fractional Caputo derivatives and Lemma 4.6 we get Developing we obtain From the fourth equation of (23) we get Substituting (33) in (32), and in view of Remark 2 and Theorem 4.1 we obtain Since we get In the other hand we have Then the maximum invariant set for is the set and according to the LaSalle’s invariance principle the DFE is globally asymptotically stable. □

Endemic equilibrium point

We investigate the endemic equilibrium point (EE) . First since does not appear in the first five equation of the system (8) then we consider and the following system For we have , so we have From (37) we getand from (38) we have Inserting (41) in (39) gives So Eqs. (35), (36) becomes Eq. (44) givessubstituting (45) in (43) we get In the other hand we have Thenand we getsubstituting (50) in (46) we obtain Suppose , then the system (34) has a unique EE denoted by , with If , then the EE denoted by of the system (34) is locally symptomatically stable when . The associated Jacobian matrix of (34) evaluated at is given by Some calculations give the following characteristic polynomialwhere If then all the eigenvalues has negatives real part. Then is locally asymptotically stable. □ We examine the global stability of the EE , for this we need the following lemma:

Lemma 3.1 [29]

Let be a continuous and derivable function. Then for any time instant If , then the EE of the system (34) is globally symptomatically stable. First all considering the following simpler model obtained by normalizing in (34) to be 1then the EE satisfy the following equations Let , and the following Lyapunov function is define for the required result Using linearity propriety of Caputo derivatives and result in Lemma 4.9, and from model (54) we have Using direct calculation, and formulas Eqs. (55), (56), (57), (58), (59), we obtain Adding up equations Eqs. (62), (63), (64), (65), we get Finally by the arithmetic–geometric means inequality, it follows thatand if in additionthen, we have . In addition we have if and only if , hence the maximum invariant set foris the singleton set and according to the LaSalle’s invariance principle the EE for the system (34) is globally asymptotically stable, whenever .  □

Numerical simulation of the model

This section of the paper presents the dynamics of the COVID-19 model (8) graphically. The main purpose is to analyze the impact of memory index and of other key parameters on the dynamics and ultimately the possible control of the current ongoing COVID-19 pandemic. The COVID-19 mathematical model (8) in Caputo sense is solved via the fractional Adams-Molten type iterative scheme. The parameters values tabulated in Table 1 which are obtained from reported infected cases in Algeria are used in the simulation results. In Fig. 3, Fig. 4, Fig. 5, Fig. 6, Fig. 7, Fig. 8 , the impact of arbitrary fractional order (i.e., memory index) is depicted graphically. The susceptible class decreases for all values of to a specific positive density as can be seen in Fig. 3. The same behavior is found for all values of . Fig. 4 depicts the dynamics of asymptomatic class for varying values of . It is observed that the peaks of the infected curves slightly decreased and occurred comparatively over longer period of time for smaller values of . The same interpretations are obtained for the population in the remaining infected classes as presented in Fig. 5, Fig. 6, Fig. 7. The dynamics of recovered or removed population for various values of is analyzed in Fig. 8. The recovered individuals curve initially increased until become stable. Furthermore, in Fig. 9, Fig. 10, Fig. 11 , we analyze the dynamical behavior of cumulative infective population () for various values of the disease transmission rates and corresponding to the population in and compartment respectively. The influence of by decreasing it with different rates to its baseline value versus the total infected individuals (i.e., ) is described graphically in Fig. 9. This interpretation is performed for four values of the memory index . A significant decrease in the peaks of infected curves is observed by decreasing the disease transmission coefficient to its estimated values as shown in Table 1. The results become more prominent for smaller values of as can be seen in Figs. 9(a-d). The impact of disease transmission rate and the fractional order versus the total infected individuals is depicted in Fig. 10(a-d). The simulation results in this case are performed by the varying values of with different rates to its baseline and for four different values of . One can interpret that the peaks of infected curves become flatten significantly with the decrease in . Finally, in Fig. 11, we perform the simulation results in order to explore the influence of the transmission rate over the total infective population. In a similar way, we depict the graphical results for distinct rates of to its baseline and for four different values of . A reasonable decrease in the infected individuals is observed with the decrease in the parameter .
Fig. 3

Dynamics of suspectable individuals for various values of .

Fig. 4

Dynamics of class for various values of .

Fig. 5

Dynamics of class for various values of .

Fig. 6

Dynamics of class for various values of .

Fig. 7

Dynamics of class for various values of .

Fig. 8

Dynamics of class for various values of .

Fig. 9

Influence of the cumulative infected individuals for (a) , (b) , (c) , (d) .

Fig. 10

Influence of the cumulative infected individuals for (a) , (b) , (c) , (d) .

Fig. 11

Influence of the cumulative infected individuals for (a) , (b) , (c) , (d) .

Dynamics of suspectable individuals for various values of . Dynamics of class for various values of . Dynamics of class for various values of . Dynamics of class for various values of . Dynamics of class for various values of . Dynamics of class for various values of . Influence of the cumulative infected individuals for (a) , (b) , (c) , (d) . Influence of the cumulative infected individuals for (a) , (b) , (c) , (d) . Influence of the cumulative infected individuals for (a) , (b) , (c) , (d) .

Conclusion

The current investigation addressed the dynamics of a new mathematical model for the transmission dynamics of the novel COVID-19 infection. To formulate the proposed model the COVID-19 infective individuals are sub-divided into detected and undetected classes. The Caputo derivative is utilized in order to better explore the disease dynamics. After the model formulation, initially, we presented the basic and necessary mathematical features of the fractional COVID-19 epidemic model. The local and global stability of both the disease-free as well as the endemic steady states are shown via considering the fractional Lyaponov functions. Moreover, with the help of non-linear least square procedure, some of the model parameters are estimated from the reported cases in Algeria while the remaining are estimated from the literature. The fractional model in Caputo sense is solved numerically via the iterative scheme of Adams-Molten type. Detailed simulation results are performed for various estimated parameters of the model. Furthermore, we also analyzed the impact of memory index (the fractional order of Caputo operator) on the disease dynamics. We believed that the current study will be helpful in curtailing the ongoing COVID-19 pandemic. The present model can be extended to more sophisticated mathematical models with a non-singular kernel in the generalized Mittag–Leffler function.

CRediT authorship contribution statement

Yacine El hadj Moussa: Conceptualization, Methodology, Resources, Investigation, Visualization, Writing - original draft, Writing - review & editing. Ahmed Boudaoui: Conceptualization, Methodology, Investigation, Visualization, Writing - original draft, Writing - review & editing. Saif Ullah: Conceptualization, Methodology, Investigation, Visualization, Writing - original draft, Supervision, Writing - review & editing, Funding acquisition. Fatma Bozkurt: Conceptualization, Methodology, Investigation, Visualization, Writing - original draft, Supervision, Writing - review & editing, Data curation, Formal analysis. Thabet Abdeljawad: Conceptualization, Methodology, Investigation, Visualization, Formal analysis, Supervision, Writing - review & editing. Manar A. Alqudah: Investigation.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
  14 in total

1.  Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission.

Authors:  P van den Driessche; James Watmough
Journal:  Math Biosci       Date:  2002 Nov-Dec       Impact factor: 2.144

2.  Modeling and simulation of the novel coronavirus in Caputo derivative.

Authors:  Muhammad Awais; Fehaid Salem Alshammari; Saif Ullah; Muhammad Altaf Khan; Saeed Islam
Journal:  Results Phys       Date:  2020-11-17       Impact factor: 4.476

3.  Asymptomatic patients as a source of COVID-19 infections: A systematic review and meta-analysis.

Authors:  Andreas Kronbichler; Daniela Kresse; Sojung Yoon; Keum Hwa Lee; Maria Effenberger; Jae Il Shin
Journal:  Int J Infect Dis       Date:  2020-06-17       Impact factor: 3.623

4.  Remdesivir for COVID-19 in Europe: will it provide value for money?

Authors:  Rafael Dal-Ré; Rita Banzi; Sophie Georgin-Lavialle; Raphaël Porcher; Reecha Sofat; Markus Zeitlinger; Frits R Rosendaal
Journal:  Lancet Respir Med       Date:  2020-12-17       Impact factor: 30.700

5.  Mathematical modeling for the outbreak of the coronavirus (COVID-19) under fractional nonlocal operator.

Authors:  Saleh S Redhwan; Mohammed S Abdo; Kamal Shah; Thabet Abdeljawad; S Dawood; Hakim A Abdo; Sadikali L Shaikh
Journal:  Results Phys       Date:  2020-11-16       Impact factor: 4.476

6.  A robust study on 2019-nCOV outbreaks through non-singular derivative.

Authors:  Muhammad Altaf Khan; Saif Ullah; Sunil Kumar
Journal:  Eur Phys J Plus       Date:  2021-02-03       Impact factor: 3.911

7.  Mathematical modeling for novel coronavirus (COVID-19) and control.

Authors:  Marei Saeed Alqarni; Metib Alghamdi; Taseer Muhammad; Ali Saleh Alshomrani; Muhammad Altaf Khan
Journal:  Numer Methods Partial Differ Equ       Date:  2020-11-26       Impact factor: 3.568

8.  A study on fractional COVID-19 disease model by using Hermite wavelets.

Authors:  Sunil Kumar; Ranbir Kumar; Shaher Momani; Samir Hadid
Journal:  Math Methods Appl Sci       Date:  2021-02-07       Impact factor: 3.007

Review 9.  The trinity of COVID-19: immunity, inflammation and intervention.

Authors:  Matthew Zirui Tay; Chek Meng Poh; Laurent Rénia; Paul A MacAry; Lisa F P Ng
Journal:  Nat Rev Immunol       Date:  2020-04-28       Impact factor: 108.555

View more
  2 in total

1.  Modeling the dynamics of coronavirus with super-spreader class: A fractal-fractional approach.

Authors:  Xiao-Ping Li; Saif Ullah; Hina Zahir; Ahmed Alshehri; Muhammad Bilal Riaz; Basem Al Alwan
Journal:  Results Phys       Date:  2022-01-15       Impact factor: 4.476

2.  Optimal control analysis of COVID-19 vaccine epidemic model: a case study.

Authors:  Arshad Alam Khan; Saif Ullah; Rohul Amin
Journal:  Eur Phys J Plus       Date:  2022-01-23       Impact factor: 3.911

  2 in total

北京卡尤迪生物科技股份有限公司 © 2022-2023.