Mathematical model of survival of fractional calculus, critics and their impact: How singular is our world?

Fractional calculus as was predicted by Leibniz to be a paradox, has nowadays evolved to become a centre of interest for many researchers from various backgrounds. As a result, multiple innovative ideas had emerged, which caused significant divisions regarding fractional calculus in the past three years. Therefore, this work is aimed at developing a mathematical model that could be used to depict the survival of fractional calculus. Six classes are herein considered to construct a mathematical model with six ordinary differential equations. All elementary analysis have been performed. Additionally, a new analysis including strength number that accounts for the accelerative information of nonlinear and linear parts of a given epidemiological model is introduced. An analysis of the second derivative of the Lyapunov function as well as an analysis of the second derivative of each class is applied to assess how a wave could be detected. It is strongly believed that this new analysis will particularly open new doors within the field of epidemiological modelling, which will aid researchers to better understand the spread of infectious diseases. The stochastic version of the suggested model was also investigated, and numerical simulations were performed. The obtained reproductive number, strength number, extinction of criticism together with numerical simulation, revealed that the field of fractional calculus will be stable will therefore have no significant effect soon.


Introduction
Fractional differentiation and integration are fast spreading topics that have become a center of interest in multiple research works due to their wider application. Researchers from all fields have since been attracted by this topic [1][2][3][4][5][6][7][8][9][10]. The topic was derived from a question raised by L'Hopital to Leibniz, which initially brought up the differentiation of exponential functions [4]. Later onwards Liouville suggested a name of the derivative with a fractional index and simultaneously warned that the properties of this new operator should not be associated with those of the classical derivative [1][2][3][4][5][6][7][8][9][10]. Thereafter, a connection of independent works of Liouville and Riemann gave birth to the well-known fractional integral, and then later to a fractional derivative. Their differential operator is a derivative of a continuous function and the power-law function t -α (1-α) . An application of the Laplace transform of this derivative leads to unusual initial conditions that would be impossible to compute [1][2][3][4][5][6][7][8][9][10]. However, these initial conditions were praised by some respected authors, simply because they can account for anomalous processes. Contrary to this, using experimental data, Caputo was able to suggest an alternative differential operator which is the convolution of the derivative of a function with power law [3]. An application of the Laplace transform to this version leads to the normal initial condition. In addition, the function was expected to be differentiable, of which some authors theoretically believed it was a restriction, even though in practice it was found to be a good differentiation, essential to evaluating changes [3]. These differential operators have been intensively studied in the last decades, causing them therefore to be both wrongly and successfully applied in so many fields. As a result, the main properties of these operators were then established. The most important property is perhaps the fact that this function exhibits a power-law process. A clear indication that such a mathematical formula cannot account for processes following fading memory crossover and many other processes like random walk and others [1,2,[11][12][13]. And a clear indication that such an operator cannot be used to model all real-world problems. Therefore, a new class of differential operators was needed, and a few variants were suggested, among which we can mention the Caputo-Fabrizio derivative and the Atangana-Baleanu derivatives with exponential and Mittag-Leffler kernels. One of the particularities of these operators is that their kernels are nonsingular [1,2,[11][12][13][14][15][16][17][18], and they do not obey some properties of classical derivative such as the index law. Thus, due to their continuous kernels, these fractional derivatives do not have similar properties like classical or fractional derivatives with the power-law kernel. This is evidence that indeed power-law, exponential, and the Mittag-Leffler functions do not play the same role in nature [10,11,19,20]. In nature, there exist problems that follow power-law processes, some that follow a decay process, and others that follow decay and then power law. In a normal and constructive field both concepts could be studied independently, because both fractional derivatives with singular kernels and fractional derivatives with nonsingular kernels are two separate subsets of fractional calculus. They should therefore be treated independently. However, instead of investigating the properties of nonsingular kernels and fractional derivatives independently, some researchers have opted to write some destructive papers mixing up these two subsets, whereby some of these papers are based on wrong analyses, birthed from the lack of understanding of the whole concept. Hence, this paper will review all issues raised against nonsingular kernels differential and integral operators. A mathematical model depicting a survival of fractional calculus based on nonsingular kernels will herein be proposed and analyzed.

Review of some anti-nonsingular kernel 2.1 Index law
One of the criticisms of differential operators is the index law, i.e., these differential operators do not satisfy the index law. Some researchers suggested that a differential operator cannot be called derivative or fractional derivative if it does not satisfy the index law. However, Atangana and Gomez argued that the index law is just a mathematical property which has many limitations to depict real-world problems with crossover behaviors [21]. Nieto and co-authors proved that there is no fractional derivative that verifies the index law [22]. Nevertheless, some researchers insisted that this law is fundamental even to replicate some real-world problems in signal analysis. Very recently, a mathematical operator was suggested in [23], and an artificial parameter was added to help satisfy the so-called index law. Nevertheless, in the paper written by Jacek and co-author, they argued that: Naturally it is not wise to conclude fractional ordinary D α t and D β t commute in their entire domain, another clear indication that fractional order derivatives will not hold even if the kernel is that of power law [23]. With the Machado-Orteguira type, we have the following parameters to be considered: the fractional order α and an asymmetry parameter θ which of course has no physical meaning. Their formula defined a derivative D α θ of the function f : R → R as Indeed the semi-group principle is naturally satisfied, namely In an attempt to apply this operator in modeling a real-world problem for which the operator was conceived, Jacek and co-authors concluded that "The generalization by Machado-Orteguira satisfied the index law principle but its phasor representation is artificial or is less natural. " [23] Here again nature has proven the limitation of some mathematical properties. The satisfaction of index law by a given operator does not always lead to good results [23]. Therefore, the argument around the index law is pointless more precisely if an artificial parameter is added as in the case of Machado-Orteguira type.

Initial conditions and zero-zero criticism
One of the arguments used was that a general Cauchy problem ⎧ ⎨ ⎩ ABC 0 D α t y(t) = f (t, y(t)), t > 0, y(0) = y 0 , t = 0, does not satisfy the initial condition if the solution exists. Indeed The geometric interpretation of the above in classical calculus is the surface under the curve of the function f (t) in [a, t] if the function f positive. The solution is of course the variation of the surface. Such an equation will have meaning only if one goes beyond the initial state which is the point t = a. At this point, the solution is the initial state or just a point nothing to do with a surface. So in general the equation ABC 0 D α t y(t) = f t, y(t) , y(0) = y 0 (5) can be viewed in the same way where The above function is continuous. Some researchers failed to understand the basic principle of nonlocality and even the concept of memory. For example, the equation The solution of the above equation is f (t) = u (t), the initial condition f (0) = u (0). Why is the initial condition not verified? The answer relies on the physical interpretation of the equation. The right-hand side f (t, y(t)) = u(t) can be viewed as a surface. Thus, the equation is the collection of information or accumulation of memory from the origin to the point t. Of course it is an accumulation that covers the area defined by u(t) = f (t, y(t)). Now, the zero-zero criticism does not hold because At t = 0, indeed we will need f (0, y(0)) = 0, but what is the physical meaning of the equation Can this equation be evaluated in t = 0 at the origin? Can we really evaluate a fractional derivative at the origin if the kernel is continuous? If yes, where is nonlocality? Of course, some authors will argue that the solution of the following equation The fact is simple at t = 0, we have f (0, y(0)) = RL 0 D α t y(t) t→0 . Within this section, we consider that the function f (t, y(t)) is defined in a negative line Let us assume that f (t, y(t)) is continuous at t = 0, then lim t→0 The right-hand sides of the equation, by replacing t with c then c → 0 + , we have Thus, if f (t, y(t)) is continuous at t = 0, the equation will also be valid at zero if f (0, y(0)) = 0, if y(t) is continuous at t = 0, if f (t, y(t)) is not continuous at t = 0. Then, indeed we have which is contradiction. This simply means at t = 0 that the equation C 0 D α t y(t) = f (t, y(t)) under the above has no solution at t = 0. There is a discontinuity at t = 0. If y(t) is not continuous at t = 0, then lim t→0 + y(t) = lim t→0 -y(t). Let us assume that f (t, y(t)) is continuous at t = 0, then lim t→0 + f (t, y(t)) = lim t→0 -f (t, y(t)) = f (0, y(0)). Again here, the equation at t = 0 will be problematic. So why is it misleading to evaluate a fractional differential equation at the origin? The answer is simple: it is a nonlocal equation that should describe memory, accumulation of past information, by evaluating it within an interval [a, a], where a is the lower boundary, the process becomes local, therefore the equation is no longer a fractional as the memory part is lost. Additionally, the equation does not have any differential part, thus cannot even be considered as a differential equation. But this can be done in the case of classical differential equations: since they are local equations, they should be evaluated at a point. Jacob and co-authors wrote a piece of paper in which they showed that equations with nonsingular kernel derivatives have a discontinuous solution at the origin. While their mathematical seems to be correct [24], the authors failed to understand a simple mathematical principle, i.e., the nonlocality. For example, they suggested that their approach [24] leads to a solution of this equation ABC 0 D α t y(t) = a, where a is constant and ABC 0 D α t is the Atangana-Baleanu derivative. Already in high school, we were taught that if a = 0, the equation may not have a solution. For example, taking a = 0, then If t = 0, then with f continuous then a = 0, which is contradicting. But if t = 0, then with t being a variable, we cannot find a solution. For example, if f (t) is a constant different from zero, then since a is a constant. It is therefore very surprising to see that a mathematician expects to find a solution for such an equation.
Here, we consider the Cauchy problem where the differential operator is the Riemann-Liouville derivative: We assume that f (t, y(t)) is continuous at t = 0. Then On the other hand, we have Integrating by parts yields Assuming that y(0) = y 0 = 0 Since f (t, y(t)) is continuous at t = 0, we have which is a contradiction. Also it is possible that Indeed there is no contradiction here because if such an equation exists then there is no solution at the origin due to the singularity. So here also the equation will not satisfy the initial condition. Whereas, if we consider the following equation: we assume that f (t, y(t)) is continuous at t = 0, then indeed The above equation shows a solution that verifies the initial condition. This is due to the fact that the above equation is a local equation thus can be evaluated at the origin.

Fundamental theorem of fractional calculus
Another point was raised against nonsingular kernel based differential operators, the fundamental theorem between the Caputo-type and the Riemann-Liouville derivative type of Atangana-Baleanu and Caputo-Fabrizio derivatives. It was reported that the nonsingular differential operators should not be used because they do not satisfy the fundamental theorem of fractional calculus. Whereas the Atangana-Baleanu in the Riemann-Liouville sense indeed is invertible with the Atangana-Baleanu integral, a clear satisfaction of the theorem. But indeed the Atangana-Baleanu derivative in the Caputo sense is not at all invertible with the Atangana-Baleanu integral as these operators are not invertible. Indeed the authors of this particular paper should not be blamed as the Caputo calculus has not been developed independently. Always, there is a great expectation that all the properties of classical differential and integration should be verified by these operators whereas both concepts are far different. However, very recently Atangana asked the question to know why there are two fractional derivatives but only one integral, whereas the two fractional derivatives were conceived differently. As an attempt to solve this problem, Atangana suggested an integral associated with Caputo type of derivatives including the Caputo integral for power-law kernel derivatives, exponential decay law, and the generalized Mittag-Leffler kernel. Indeed they all verified the fundamental theorem of fractional calculus with their respective derivatives. In general, the Caputo integral of a continuous function is given by where RL 0 J α t is the Riemann-Liouville integral. The author of [12] suggested two types of fractional calculus including Caputo calculus and the Riemann-Liouville calculus. Indeed, criticisms are good when they aim at developing or enhancing science; however, when they are not based on truth or have different aims, they can be classified as virus as they can cause harm within a given field. Therefore, in the following section, we will introduce a mathematical model about the spread of harmful criticism within the field of fractional calculus.

Mathematical model of anti-nonsingular kernel
In this section, we construct a mathematical model that could depict the effect and impact of harmful criticism within the field of fractional calculus. Let be the in-flow, the number of individuals that are recruited in the field of fractional calculus. F c (t) is the class of researchers using fractional differential operators with nonsingular kernels susceptible to be mislead by the harmful criticisms. I(t) is the class of researchers that have been affected by the harmful papers. I P (t) is the class of researchers that have been affected but still have positive opinions about nonsingular kernel derivatives. I N (t) is the class of researchers that have been affected and have negative opinions about nonsingular kernel derivatives. R(t) is the class of researchers that overcome divisions criticism. D(t) is the class of researchers that die or retire or leave fractional calculus due to division. The field is very attractive as the concept is widely applicable, therefore, we will assume that a new researcher will join the field, thus in this model, we consider as the recruitment rate. Naturally, some researchers working within this field may die, thus they will not be working within the field, we will then consider d as the removal rate. Indeed d can also account for the retirement. β will be considered as a coefficient of contact, a proportion that a re- searcher reads a paper with a negative content about nonsingular kernel derivatives. φ 2 is the rate at which a researcher working in fractional calculus joins class D due to division. τ is the rate at which who are affected by criticism, but still have positive opinions about nonsingular kernel derivatives. δ is the recovery rate. ψ is the rate with which a researcher joins I N class, γ 1 is the rate at which individuals with positive opinion join class D due to division. φ 1 is the recovery rate of I P class. σ is the rate at which class I N joins class D due to division. γ 3 is the recovery rate of class I. γ 2 is the recovery rate of class D. In the Fig. 1, we present a mathematical diagram depicting the dynamic described.
In this model, we have assumed that no recovered person joins the class F c just to see if the class F c will survive even without researchers that were impacted by criticisms. A mathematical model associated with the above is given as follows: with the initial conditions

Positiveness and boundedness of solutions
In this subsection, we present a detailed analysis underpinning the conditions under which positiveness the solutions of the suggested model for different types of differential operators including integer and non-integer orders holds. To achieve this, we write For the function I P (t), one can write ≥ -(d + γ 1 + φ 1 )I P , ∀t ≥ 0, since I(t) ≥ 0, ∀t ≥ 0. From above, we get Doing the same routine for other classes, the following inequalities are obtained: We shall define the norm where D g is the domain of g. Using the above definition, the following inequality can be achieved for the function F c (t): This yields

Positive solutions with nonlocal operators
In this section, we prove the positiveness of solutions for a fractional calculus model with nonlocal operators. If all the initial conditions are positive with nonlocal operators, all solutions are positive. To do this, we start with the Caputo-Fabrizio case , ∀t ≥ 0, , ∀t ≥ 0, , ∀t ≥ 0, , ∀t ≥ 0, With the Caputo derivative, we obtain With the Atangana-Baleanu derivative, we get , ∀t ≥ 0, , ∀t ≥ 0, , ∀t ≥ 0, , ∀t ≥ 0, If we do the same routine for fractal-fractional operators [17], for a power-law kernel we can have the following: here c is a time component. With an exponential kernel, we can write , ∀t ≥ 0, , ∀t ≥ 0, , ∀t ≥ 0, , ∀t ≥ 0, With a Mittag-Leffler kernel, the following can be done: , ∀t ≥ 0, , ∀t ≥ 0, , ∀t ≥ 0, , ∀t ≥ 0,

Analysis of equilibrium points
In this subsection, a detailed analysis of equilibrium points is presented. We start with the disease-free equilibrium To obtain the equilibrium points, we solve the following system: , which are disease-free equilibrium points. We should not consider that case, we therefore assume I * = 0 such that , ,

Reproduction number
To obtain the reproduction number, we consider the following equations: Now, we have to calculate the matrices F and V by employing the next generation matrix approach as follows: Here, Indeed, we evaluate at disease-free equilibrium to have In this case, we can have the following: Then, the reproductive number can be By evaluating F at E * , we have Then, we can evaluate the following parameter: Before we present a detailed derivation on stability analysis of the suggested model, we need first to evaluate the Jacobian matrix of the system at the equilibrium points. We first where l 4 = d +γ 2 . The associated eigenvalues to J are -d, -d, -l 1 , -l 2 , -l 3 , -l 4 and βd-dl 1 -l 1 φ 2 d+φ 2 . We can see that all eigenvalues have a negative real part except βd-dl 1 -l 1 φ 2 d+φ 2 . Nevertheless, rearranging this last one, we get The above expression can be negative if βd l 1 (d+φ 2 ) < 1. This leads to In the following figures, we present numerical simulations of the reproductive number for different values of theoretical parameters including beta and delta. The numerical simulations are depicted in Figs. 2, 3, 4, 5, 6, and 7. The numerical results confirm that the reproductive number of such a model is always less than 1, which is a clear indication of less impact, perhaps due to the fact that published papers against nonsingular kernels are not well-grounded either theoretically or in practice. This can be also justified with the fact that thousands of papers using derivatives with nonsingular kernels have been published with outstanding results in both theory and applications.

Strength number
In the last decades, the concept of reproduction has been employed intensively in epidemiological modeling as it has been recognized as a useful mathematical formula to evaluate reproduction in some given infections disease. As the theory suggested, one will find two component F and V , then will be used to reproduce the reproductive number [25]. The component F is very interesting as it is obtained by deriving the nonlinear part of the infected classes.

∂ ∂I
At the disease-free equilibrium point, we have (62) In this case, we can have the following: Then A 0 = 0 means the spread will not have a renewal process, therefore will have a single magnitude and die out. A 0 > 0 means there is as strength that will lead to renewal process, therefore the spread will have more than one wave. However, a clear interpretation of the above number will be provided by biologists.

Stability analysis
To discuss the stability of the model, we need first to evaluate the Jacobian matrix of the system at the equilibrium points of the model. First, we study the researcher free criticism equilibrium point and its stability The eigenvalues associated with J are -d,-Q 2 ,-Q 3 ,-Q 4 ,-dφ 2 , and βd-dQ 1 We can see that all the eigenvalues have negative real parts except βd-dQ 1 -Q 1 φ 2 d+φ 2 . If we rearrange the last eigenvalue as λ 6  Proof To prove theorem, we consider the Lyapunov function as follows:

Lyapunov for endemic case
Therefore, applying the derivative with respect to t on both sides, we get the following: Now, we put values in the above equation for derivatives we can have the following: The above equality can be rearranged as For simplicity, we can rewrite the above equality as follows: It can be easily seen that if < , this yields dL dt < 0, however We follow that the largest compact invariant set for the suggested model model in is the point {E * }, the endemic equilibrium of the considered model. From Lasalle's invariance concept, we can conclude that E * is globally asymptotically stable in if < . The above derivation does not sometimes provide a clear condition on the sign of the first derivative of the Lyapunov function despite being used in many different papers. While the derivation is very useful, more efforts are required to indeed access the sign of the first derivative of the Lyapunov function. However, there is an alternative way to obtain a clear condition using the equilibrium points, and it will be presented below to show readers the difference between the two approaches. The following are the results for the model at steady state: then endemic equilibrium is globally asymptotically stable.
Proof Consider the following nonlinear Lyapunov function for the suggested model, so that the corresponding unique endemic equilibrium point E * 1 exists by letting R 0 > 1: The time derivative of above is Using the solution of the system gives Putting all in the above after simplification gives Since the arithmetic mean exceeds the geometric mean, we have the following interpretation: Further, if the following inequality holds then L (t) ≤ 0 for R 0 . Thus, L(t) is a Lyapunov function in . Therefore, by Lasalle's invariance principle, we have Thus, every solution of the model tends to its unique endemic equilibrium for the associated reproduction number as t → ∞. The endemic equilibrium point E * 1 is global asymptotically stable whenever R 0 > 1.

Lyapunov function for the disease-free case
Theorem 3 For the survival of fractional calculus model, the harmless equilibrium or the critical equilibrium is globally asymptotically stable within the stable feasible if R 0 < 1 and unstable if R 0 > 1.
Proof We use the following: Taking its derivative with respect to t, we have Putting all from the system, we get · L = βF c I Nl 1 - ≤ (R 0 -1)(I + I P + I N ).
Therefore if R 0 -1 < 0, then Lyapunov decreases · L(t) < 0 since we have shown that I + I P + I N is positive. Therefore L is a Lyapunov function within the feasible biological interval, and the bigger compact invariant set in {I, I P , I N ∈ : dL dt ≤ 0} is the point E 0 . According to Lasalle's invariance concept [26], each solution with initial condition in leads to E 0 when t → ∞ only if R 0 ≤ 1. Thus, one can conclude that the disease-free equilibrium E 0 of the fractional calculus model is globally asymptotically stable.

Second derivative of Lyapunov
The global stability of endemic equilibrium points is evaluated by means of the first derivative of the Lyapunov function. Without loss of generality, the analysis of the first derivative provides us with useful information that can be completed by the analysis of the second derivative. For example, the second derivative provides us with the curvature accordingly to its sign, while the first derivative of these Lyapunov functions gives us some useful information about the spread of the disease. We strongly believe that its second derivative could possibly add more information. Here, Then we have and · R(t), and · D(t) with their respective formula and then putting together positive and negative factors, we have Then If 1 > 2 then d 2 L dt 2 > 0, Then the interpretation associated with the sign of second order follows.

Equilibrium points for second order
As discussed before, the second derivative provides very useful information on curvatures.
Also the equilibrium points of such can also provide very useful information. We shall now present the equilibrium points of the second order derivative of our solutions.
At the disease-free equilibrium, we have the following solution For endemic equilibrium, we obtain Although equilibrium points have been used successfully to give some information about the spread of any disease, it is worth noting that these points are obtained using only the first derivative, which of course just gives an indication about the rate of change. Whereas F * * c , I * * , I * * P , I * * N , R * * , and D * * are able to give some information about infection points, which can be used to detect wave. Therefore, if indeed we have two equilibrium points, then one could think about two waves. Indeed, on any graph, the second derivative will correspond to the curvature, or in other words concavity, of the graph.
We now evaluate the sign of I * * , I * * P , and I * * N . As prescribed before, we will check the conditions under which we can have I * * > 0, I * * = 0, and I * * < 0. However, where by definition Dividing by I N , owing to the fact that the class I(t) is not zero everywhere, yields By definition F c N < 1 also I N < 1, we have and we find Therefore, if the above criteria are satisfied, ·· I < 0 the function has a local maximum otherwise then ·· I > 0, the function has a local minimum or ·· I = 0, there is an infection point that can be found. On the other hand, we write Since by definition I > I P , then τβIF c N τ I(l 1 + l 2 ) + l 2 2 I P < 0, τβIF c N + l 2 2 I P < τ I(l 1 + l 2 ), τβIF c N + l 2 2 I < τ I(l 1 + l 2 ) + l 2 2 I, Also, we write Therefore, to have a maximum value on time for the three classes, we need

Existence and uniqueness
In this section, we present a detailed analysis about the existence and uniqueness of the system of equations describing the survival of fractional calculus. To achieve this, the following theorem is to be verified.

Theorem 4 Assume that there exist positive constants
We now recall our model We start with the function F 1 (t, F c , I, I P , I N , R, D). Then, we will show that Then, we write where κ 1 = {2d 2 + 2β 2 I 2 ∞ }.
≤ κ 2 I 1 -I 2 2 , We verified the first condition for all function. We now verify the second condition for our model.
∞ under the condition d 2 L < 1. Then, the solution of our system exists and is unique under the condition

Stochastic model of survival of fractional calculus
Since there is no clear reason as why the division was initiated within the field of fractional calculus, we assume some environment noises that could be characterized by envy, lack of understanding of the concept of fractional derivatives with nonsingular kernels, and different human behaviors that could lead to divisions. The suggested mathematical model will be converted to a stochastic model as follows:

Existence of unique global positive solution
In this subsection, we present the existence of a unique positive solution of the suggested model.  D(t)) of the stochastic model on t ≥ 0, and the problem solution will maintain in R 6 + with unit probability.
Proof As the coefficient of the equation is locally continuous in the Lipschitz sense for the given initial size of population (F c (0), I(0), I P (0), I N (0), R(0), D(0)) ∈ R 6 + , there must exist a unique solution (i.e., local solution) (F c (t), I(t), I P (t), I N (t), R(t), D(t)) on t ∈ [0, τ e ), where τ e denotes the explosion time. In order to show that actually the solution is global, one has to prove that in fact a.s. τ e = ∞. Let us consider a positive real number k 0 and large enough so that all of the initial values of the states lie within { 1 k 0 , k 0 }. Further, let us define the stopping time for each nonnegative integer k greater than or equal to k 0 . We assumed here that inf φ = ∞ whenever φ denotes the empty set. By looking into the definition of stopping time, one can say that τ k is monotonically increasing k → ∞. Set If, for all nonnegative values of t, we show that τ ∞ = ∞ a.s., then we can say that τ e = ∞ and a.s. (F c (t), I(t), I P (t), I N (t), R(t), D(t)) ∈ R 6 + . Thus, we have to prove that τ e = ∞ a.s. If the conclusion is assumed to be false, then there must exist two constants 0 < T and ∈ (0, 1) such that Next, we will define a function H : R 6 + → R + from the C 2 space such that H(F c , I, I P , I N , R, D) = F c + I + I P + I N + R + D -6 (135) -(log F c + log I + log I P + log I N + log R + log D).
In above, H : R 6 + → R + may be defined through the following relation: and LH(F c , I, I P , I N , R, D) Here, the formulation of K shows that it is positive and independent state variable as well as independent variable. Therefore Integrating both sides of the above equation from 0 to τ k ∧ T, we have Setting k = {T ≥ τ k } for k 1 ≤ k and thus P( k ) ≥ . Note that, for each w in k , there must exist at least one F c (τ k , w), I(τ k , w), I P (τ k , w), I N (τ k , w), R(τ k , w), D(τ k , w), which equals 1 k or k. Hence (F c (τ k ), I(τ k ), I P (τ k ), I N (τ k ), R(τ k ), D(τ k )) is not less than k -log k -1 or log k -1 + 1 k . As a result, Here, the notation 1 w represents the indicator function of . Letting k → ∞ will lead to the contradiction ∞ > H(F c (0), I(0), I P (0), I N (0), R(0), D(0)) + TK = ∞, which implies that τ ∞ = ∞ a.s. and completes the proof.

Extinction for criticism
In this subsection, we aim at deriving the conditions for which criticism could be extinct.
In order to avoid incomprehension by some readers, we define some well-known notations [27]. Let us consider We also define the relation depending on the threshold R 0 for the proposed model of survival of calculus . (146)

Lemma 1
We assume that the classes (F c , I, I P , I N , R, D) represent a system of the solution of our model with initial data The methodology used to establish the above lemma can be found in Zhao and Jiang's work, readers are welcome to check references [27]. To avoid repetition, we will omit the proof here as it is similar to what has been done in [27][28][29][30][31][32].

Theorem 6 Under the condition that d >
That is, I(t) → 0 exponentially, which implies that the process will die out and the concept of nonsingular derivative will survive as well as fractional calculus with unit probability. Additionally, .
Proof Applying integral on both sides and dividing by t, we obtain Now, applying Ito's formula to the above system, we obtain Integrating the above within the interval [0, t] and dividing the same by t yields Additionally, which is local continuous martingale and M(0) = 0. By the first lemma and t → ∞, we obtain If R 0 < 1 is satisfied, then The above inequality implies that lim t→∞ I(t) = 0. Also, by applying the same principle on I P (t), we have Similarly, Also the above yields In a similar way, we can show that However, to achieve the value of lim t→∞ F c (t) , we apply integral on both sides and divide by t to obtain Thus, taking the limits when t → ∞, we get Similarly, we have Thus, applying the limit as t → ∞, which completes the proof. We now present the extinction of the classes I(t), I P (t), and I N (t) when the fractional derivative is the Caputo-Fabrizio version By dividing with t, the above produces Since F c (t) represents the class of researchers working in fractional calculus, clearly, for Considering the fact that F c N ≤ 1 and that F c is bounded, we get For the I(t) class, we have the following: and dividing by t yields We write Thus, we have Similarly, we write and dividing by t yields We write However, lim t→∞ I(t) = 0, thus lim t→∞ I P (t) = 0.
Applying a similar procedure on I N and R also, we obtain However, for D class, we have . .
We now present a discussion underpinning the extinction of species where the model is with fractional derivative in the Caputo sense.

Lemma 2 Assuming that a function x(t) is continuous and bounded for t ∈ [0, ∞), then for
For proof, since x(t) is continuous and bounded, then We now assume that x(t) is bounded, then We now discuss the condition of existence of the following operator for any 0 < α ≤ 1. To do this, we shall present the following definition: Then Dividing by t α , we get Then, we have We consider the class I, we have the following: and dividing by t α yields Since ∀t ∈ [0, T] F c < N , then and Thus we write Then, we get Indeed l 1 + σ 2 2 2 > 0, also 1 -R 0 > 0 since it was assumed that R 0 -1 < 0, then Using the same routine, we have We write However, lim t→∞ C I α = 0, therefore In the same way, Finally, we evaluate the class D Thus we write Therefore, we write With the above, we can now show the extinction with ABC derivative Dividing by t α yields and The left-hand side gives zero, then For the class I(t), we have Dividing by t α yields

Optimal control for the model of survival of fractional calculus
In this section, we present optimality conditions for the model of survival of fractional calculus by using Pontryagin's maximum principle [33]. To achieve this, we will rearrange the considered model by using four control variables as four possible control strategies. The control variable u 1 is the discussion about singular and nonsingular kernels in social media such as Researchgate, Facebook. The control u 2 describes the number of books, papers accepted about singular and nonsingular kernels. The control u 3 describes the fairness of editorial and publishers. The control u 4 is a conference with good talk where fractional calculus is discussed. By addition of the control functions mentioned above, we modify our model as follows: The objective functional can be represented by The parameters q 1 , q 2 , q 3 , q 4 , a 1 , a 2 , a 3 , a 4 are the weighted parameters. The existence of the control functions [15] can be ensured by the assumptions the set of U is nonempty, convex, bounded, and closed, the Lipschitz property of the state system, and the convexity of the integrand of the objective functional with respect to the controls on the set U. Based on Pontryagin's maximum principle, we establish the Hamiltonian H given by H = a 1 u 2 1 + a 2 u 2 2 + a 3 u 2 3 + a 4 u 2 4 + q 1 F c + q 2 I + q 3 I N + q 4 I Pq 5 R Then, we obtain the following necessary conditions: , with the transversality conditions λ k (t f ) = 0 for k = 1, 2, 3, 4, 5, 6, and control variables are given by Thus, the optimality conditions can be achieved as u * 4 = min u 4 , max 0, 9 Numerical solution of the model with exponential decay process In this section, we present a numerical scheme based on Newton polynomial to solve the suggested mathematical model for different differential operators [16]. We start with the Caputo-Fabrizio case for a numerical solution of the mathematical model of survival of fractional calculus with nonsingular kernel where X = (F c , I, I P , Applying the associated integral, interpolating f i (t, X i ) using the Newton polynomial within [t n , t n+1 ] yields 23 12 F(t n , X n ) t - 4 3 F(t n-1 , X n--1 ) t + 5 12 F(t n-2 , X n-2 ) t

Numerical solution of the model with Mittag-Leffler process
In this section, we present a numerical solution for the model where the time derivative is the fractional derivative with the generalized Mittag-Leffler kernel. Here again, we adopt the numerical scheme based on the Newton polynomial interpolation [16]. Therefore, applying the scheme on the model, we obtain the following algorithm. It is important to note that a model with the generalized Mittag-Leffler helps capture processes with a passage from stretched exponential to power law with no steady state: where X = (F c , I, I P , Applying the associated integral, interpolating f i (t, X i ) using the Newton polynomial within [t n , t n+1 ] yields (nj + 1) α 2(nj) 2 + (3α + 10)(nj) +2α 2 + 9α + 12 -(nj) α 2(nj) 2 + (5α + 10)(nj) +6α 2 + 18α + 12

Numerical solution of the model with power-law process
It has been presented in the literature that Caputo derivative, which is based on the powerlaw kernel, is suitable for modeling real world problems exhibiting power-law processes. Therefore in order to include in our model the effect of power law that could be followed by the dynamic observed in the community of fractional calculus, we replace the time derivative with the Caputo derivative. Again here, we adopt a numerical scheme based on the Newton polynomial interpolation [16]: where X = (F c , I, I P , Applying the associated integral, interpolating f i (t, X i ) using the Newton polynomial within [t n , t n+1 ] yields F t j-1 , X j-1 -F t j-2 , X j-2 + α( t) α 2 (α + 3) n j=2 F t j , X j -2F t j-1 , X j-1 + F t j-2 , X j-2 .

Numerical solution of the model with fractal-fractional exponential kernel
Fractal-fractional differential operators have been suggested to depict processes following fading memory with some self-similar patterns. In this section, in order to include into the mathematical model the effect of fading memory with self-similar patterns, we replace the time fractal-fractional derivative with the exponential-decay kernel. We adopt a numerical scheme based on the Newton polynomial interpolation [16]: where X = (F c , I, I P , I N , R, D), F(t, X) = f i (t, X i ) i=1,...,6 .

Numerical solution of the model with fractal-fractional Mittag-Leffler kernel
Processes exhibiting a passage from stretched exponential to power law with no steady state and self-similarities could be modeled using a fractal-fractional derivative with the generalized Mittag-Leffler kernel. In this section, in order to make our model more complex, we present the numerical solution when the time derivative is replaced with the fractal-fractional derivative with the generalized Mittag-Leffler kernel. A numerical scheme based on the Newton polynomial interpolation [16] is used: Applying the associated integral, interpolating f i (t, X i ) using the Newton polynomial within [t n , t n+1 ] yields

Numerical solution of the model with fractal-fractional power-law kernel
Real world problems exhibiting power-law processes with self-similar patterns could be modeled using a differential operator based on the fractal-fractional with power-law kernel. In this section, we replace the time derivative with the power-law fractal-fractional derivative. The numerical solutions are obtained using a numerical scheme based on the Newton polynomial interpolation [16]: where X = (F c , I, I P , I N , R, D), F(t, X) = f i (t, X i ) i=1,...,6 . (243) Applying the associated integral, interpolating f i (t, X i ) using the Newton polynomial within [t n , t n+1 ] yields X n+1 = ( t) α (α + 1) The class of researchers affected by criticisms Figure 10 The class of researchers affected but still with positive opinion AB 0 D α t R(t) = (φ 1 I P -dR + δI + γ 3 I Nγ 2 R) + σ 5 G 5 (t, R)W 5 (t), AB 0 D α t D(t) = (γ 2 R + γ 1 I P + φ 2 F c + dI N -dD) + σ 6 G 6 (t, D)W 6 (t) with the initial conditions (246) Figure 11 The class of researchers affected but with negative opinion (247) Figure 13 The class of researchers that die, retire, or leave fractional calculus Figure 14 The class of researchers susceptible for σ 1 = 0.000011 The obtained results from deterministic and stochastic models are in perfect agreement with the obtained reproductive and strength numbers.

Conclusion and prediction
In the last two years, concerns have been raised by publishers, editors in chiefs, and younger researchers regarding the division created within the field of fractional calculus. The division occurred due to the introduction of fractional differential operators with nonsingular kernels, which indeed brought a revolution within this field, since these new dif- Figure 15 The class of researchers affected by criticisms for σ 2 = 0.00031 Figure 16 The class of researchers affected but still with positive opinion for σ 3 = 0.00032 ferential and integral operators opened new doors both in theory and application of fractional calculus to solve real life problems. However, for some reasons, some researchers expected differential operators with nonsingular kernels to satisfy the properties of classical derivative just as fractional derivative with singular kernel, which is absolutely impossible because these operators are not used for the same purposes. For this reason, a new mathematical model depicting the survival of fractional calculus field was developed. New analysis that could be more informative, for example, toward the spread of any in- Figure 17 The class of researchers affected but with negative opinion for σ 4 = 0.00033 Figure 18 The class of researchers that overcome criticisms for σ 5 = 0.00067 fectious disease, was herein introduced. This analysis included the strength number that is obtained by taking the second derivative of a nonlinear part of infectious classes and then applying similar routine to obtain the reproductive number. The utilization of second derivative of the Lyapunov function was introduced to find the sign of the second derivative of each class, with the aim to detect waves. The mathematical model was further converted to a stochastic one and solved numerically using the numerical scheme based on the Newton polynomial interpolation. The results obtained from determinis-