A compartmental fractional-order mobbing model and the determination of its parameters

In this study, a mathematical model is presented with fractional-order differential equations in the Caputo sense that examines the time-dependent changes of variables representing individuals who are exposed to mobbing, individuals who are not exposed to mobbing, or individuals who gained resistance to mobbing, and individuals who practice mobbing. Existence and uniqueness, and bound-edness and non-negativity of the solutions of the proposed model are examined. Additionally, the data set containing the time-dependent changes of these variables has been used as a basis to examine the dynamics in a population. By the data set, the approximate results of 9 different parameters used in the proposed mathematical model with ODE have been found with the lsqcurvefit function. The parameters obtained from the ODE system are rewritten into the proposed fractional model, and the value of the derivative order that gives the minimum error has been investigated. The related Matlab codes are presented in the study. Also, while it is observed that a decrease occurs in the number of individuals exposed to mobbing, there is an increase in the number of individuals who are not exposed to mobbing or who gain resistance to mobbing and individuals who practice mobbing. All the obtained results are shown in graphical detail.


Introduction
The fractional derivative and the fractional integral are the more general forms of the known integer derivative and integral [1].The history of the fractional derivative and integral is very close to the history of the integer derivative and integral.However, when we look at recent history, these have not been given the importance they deserve due to most mathematicians being unfamiliar with fractional order [2].The basic advantages of fractional derivatives are flexibility and non-locality.Nonlocality explains that the fractional order derivative has a memory effect, which means that future states depend on the present state as well as on past states [3].Because these derivatives are of fractional order, they can approximate real data with more flexibility than classical derivatives [4,5].In 1695, many mathematicians such as L'Hospital and later Lacroix, Liouville, Laurent, Riemann, Caputo, and Atangana worked in this field and each made their own definitions of fractional derivative.The most well-known and used of these definitions are Caputo and Riemann-Liouville definitions.The Caputo derivative has some advantages over the Riemann-Liouville derivative [6].Some of these can be summarized as follows: The Caputo derivative used in this study is very useful when dealing with a real-world problem because it allows conventional initial and boundary conditions to be included in the formulation of the problem, and in addition the derivative of a constant is zero [7].Fractional-order differential equations (FODEs), being a form of interest to scientists in mathematical modeling in recent years, are used in the fields such as physics [8], thermodynamics [9], viscoelastic [10], electrical theory [11], mechatronics [12], medicine [13], chemistry [14,15], chaos theory [16], finance and economics [17], biological systems [18,19], etc.The term mobbing was first used by Leymann (1990) to describe abusive behavior in the workplace [20].Workplace mobbing is a phenomenon that has actually only been identified since the 1990s and has received scientific and legal attention, especially in Europe, but has recently become a subject of increased attention in the United States and Canada.Leymann determined two criteria for an action to be considered mobbing: These are; The action must occur very frequently, at least once a week, and must continue for a long period of time, at least 6 months [21].Mobbing, which is frequently encountered in workplaces, is a situation that has existed since the beginning of business life and which can be called as the pressure exerted by other individuals or superiors in the business environment to expel an individual who is seen in the business environment [22] As it can be seen in all business environments, regardless of language, race or belief, it is a serious issue that every individual, regardless of gender, can face.Mobbing seen in business life can be seen in different types such as blackmail, pressure, and intimidation.Mobbing affects the individual's mental health and sense of well-being and negatively affects the organization in which it occurs [23].Minimizing and resolving these conflicts is of vital importance for business managers in order to achieve the targeted gains in the workplace.Because these disagreements between workplace employees will cause a decrease in working efficiency [24].Over time, jealousy in the working environment and the effort to be better can bring many problems.So much so that due to these and similar reasons, it is inevitable that a person's colleague or management puts pressure on the employee [25].Today, mobbing is considered a violent crime in many European countries, especially Belgium, Sweden, and Finland.In addition, in the USA, psychological harassment is considered within the scope of 'workplace violence' along with physical harassment in many state laws, but in practice, it is punished with mobbing only when it is related to issues such as gender, race, and religious discrimination.Workplace violence has become an important issue in Turkey in recent years too.However, there are limited studies on mobbing in Turkey.Moreover, there are no statistics on the different forms of mobbing in the workplace, nor is there any official recognition of its existence.Despite these, some statistics have been obtained in recent years.For example, 31.1% of healthcare workers in Turkey have encountered mobbing in the last year, and the mobbing frequency of 48.6% of these workers is 1 to 3 times a year [26,27].Additionally, between 2011 and 2016, 38262 people called the Ministry of Labor and Social Security's line number 170 to get information, ask for psychological support, and mobbing situations.It was reported that 57% (21,922) of these calls came from male employees and 43% (16,340) from female employees [28].
Groups of individuals with certain characteristics in a closed community always interact with each other in the form of compartments.Such interactions can be modeled in mathematics with a structure based on the SIR model [5].Thanks to this modeling, future predictions can be made and control and prevention strategies can be developed.There are many mathematical modeling studies in the literature that examine interactions through compartments by using FODEs.For some of these, see [29,30].Employees who are exposed to mobbing in a workplace always interact with other employees and employees who practice mobbing.In this study, the time-dependent changes of these compartments were modeled with the help of fractional differential equations.The rest of the article continues as follows.
• Section 2 includes some basic definitions and concepts used in this study.
• The fractional-order mobbing mathematical model is presented in Section 3, and a qualitative analysis of the model is performed.

Preliminaries and definitions
In this section, the basic definitions and theorems used in defining the proposed mathematical model and in its qualitative analysis are included.
We have adopted a dynamic system defined by the Caputo fractional operator.The system is where . ., n is satisfied.Let us consider J as the Jacobian matrix of system (2).If the eigenvalues λ i for i = 1, 2, . . ., n obtained from the equation satisfy Routh-Hurwitz stability conditions [32,33] or inequalities given as then the equilibrium point f is locally asymptotically stable (LAS) for system (2) [7].

Fractional-order mobbing mathematical model
The dependent variables used in the model consisting of fractional-order differential equations in the meaning of Caputo to be proposed, where t represents the time-parameter, are given in Table 1.The number of individuals who are sensitive to mobbing in a workplace at time t Y(t) The number of individuals who were not exposed to mobbing or became resistant to mobbing in a workplace at time t Z(t) The number of individuals applying mobbing in a workplace at time t In addition, the parameters used in the model and whose approximate values will be determined in numerical studies are given in Table 2.The rate of employees who were exposed to mobbing but became resistant to mobbing after a while σ 1 The rate of leaving the workplace of individuals who are exposed to mobbing as a result of the interaction of individuals who are exposed to mobbing with individuals who are not exposed to mobbing or gained resistance to mobbing σ 2 The rate of leaving the workplace of individuals who are not exposed to mobbing or gained resistance to mobbing as a result of the interaction of individuals who are exposed to mobbing with individuals who are not exposed to mobbing or gained resistance to mobbing κ The increase rate of individuals who apply mobbing η The carrying capacity of individuals who apply mobbing The proposed mathematical model consisting of FODEs in the meaning of Caputo is where The parameters in (5) satisfy the following inequality: In addition, system ( 5) is defined with positive initial conditions X(t 0 ) = X 0 , Y(t 0 ) = Y 0 , and In this sense, the model has the following assumptions: • Every individual who is new to work has the potential to be subjected to mobbing by individuals who practice mobbing.• Individuals, who practice mobbing, increase according to the logistic growth rule.
• Individuals who are exposed to mobbing and individuals who are not exposed to mobbing or have become resistant to mobbing, interact with each other.• Employees exposed to mobbing either leave their jobs or continue working.These employees either gain resistance to mobbing or over time they no longer care about mobbing.• Insufficiency of factors such as salary and working environment does not affect individuals who practice mobbing.
In terms of the scenario of the model, there is a ratio in the total population of the population size of each compartment.Moreover, since these ratios (< 1) are numerically small enough, it will make it easier for the optimization metric value to be smaller than the options value in terms of parameter analysis in numerical studies.It is also obvious that the ratios multiplied by the total population size will be the population size of each compartment.Let are applied to system (5), then we have where As a result of transforms in (7), x(t), y(t) and z(t) represent the proportion of each compartment within the total population.The above-mentioned scenario of the interactions of the variables regarding the parameters used in system ( 8) is shown in Figure 1.

Existence and uniqueness of the solutions
The existence and uniqueness of the solutions of Eq. ( 8) is discussed in this section.
Proposition 1 There exists a unique solution of fractional-order system (FOS) in Eq. ( 8) with each non-negative initial condition.
Proof Existence and uniqueness of system (8) will be displayed in the region Ω×(0, T] such that For proof, it is followed the manner presented in [34]. Considering any X and X, the followings are obtained from (11).Therefore, we have where Hence, it is found ∥G(X) − G(X)∥ ≤ L∥X − X∥, such that L = max(φ 1 , φ 2 , φ 3 ).G(X) satisfies the Lipschitz condition.Therefore, it is shown the existence and uniqueness of the solutions of Eq. ( 8).

Boundedness and non-negativity
In this section, the solutions of Eqs. ( 8) are non-negative and boundedness because they have interacting population densities.

Proposition 2
The solutions of Eqs.(8) starting in R 3 + are uniformly boundedness and non-negative.
Proof For the proof, the dealing presented in [35] was used.From the third equation of the system (8), we have If κz is added to both sides of this inequality, then it is and so, By using the standard comparison theorem for FODEs, where E α is the Mittag-Leffler function.Considering Lemma 5 and Corollary 6 in [36], we have z(t)≤η, t→∞.Let us consider as x(t) + y(t) = M(t).Therefore, we have and so, In this case, one has M(t)=M(0)E α (−β(t) α )+π(t) α E α,α+1 (−β(t) α ).Therefore, we have M(t)= π β for t→∞.As a result, it is convinced that all solutions are confined to the region Ω, in which Consequently, the solutions of FOS starting in R 3  + are uniformly bounded in the region Ω.Lastly, the solutions of Eqs. ( 8) are non-negative is analyzed.From (8), the followings are found.With respect to Lemmas 5 and 6 in [37], these solutions are non-negative.■

Equilibrium points
The equilibrium points E(x, y, z) of system ( 8) are obtained by solving the equations given as Therefore, we have the system From the third equation of system (12), z= 0 or z=η is found.
i. Firstly, let z= 0. From (12), it is If the equilibrium component y = γx (β+σ 2 x) obtained from the second equation of ( 13) is substituted in its the first equation, the polynomial, is obtained.Because of (10), it is clear that the coefficient of x 2 is positive, and the constant term is negative.In this case, the discriminant of ( 14) is positive, and so, it has two roots with opposite signs.Therefore, it is obtained the equilibrium point E 0 (x, γx (β+σ 2 x) , 0), always existing where The components of the equilibrium point E 0 are not negative, and this point always exists positively due to (10).
ii.Let z * =η for the positive equilibrium point represented by E 1 (x * , y * , z * ).In this case, the system, is found by Eqs.(12).From the second equation of ( 16), y * = (νη+γ)x * β+σ 2 x * as the equilibrium component is found.Substituting this component in the first equation of ( 16), we have the polynomial where Because of (10), it is clear that a 1 > 0 and a 3 < 0. Therefore, regardless of the sign of a 2 , the discriminant of ( 17) is positive and its roots are real.Since the multiplication of roots is negative due to a 3 a 1 < 0, the positive equilibrium point is E 1 (x * , (νη+γ)x * β+σ 2 x * , η) where It is clear that the equilibrium point E 1 exists in any case due to (10).

Stability analysis
The Jacobian matrix obtained from system (8) by partial derivatives is i.By calculating the jacobian matrix at the equilibrium point E 0 (x, γx β+σ 2 x , 0), we have If the eigenvalues found from the equation Det(λI 3x3 − J(E 0 )) = 0 are calculated, one of the eigenvalues is found as λ 3 = κ.It is clear that λ 3 > 0 due to Eq. (10).Therefore, the equilibrium point E 0 is an unstable point according to the Routh-Hurwitz stability criterion and inequality (4) without the need to calculate other eigenvalues.
Consequently, the equilibrium points can be accessed in Table 3.

Numerical studies
In this section, the qualitative analysis results obtained from our proposed model are illustrated by numerical simulations.For college A, both the most ideal derivative order for the model and the most suitable parameters are shown in inequality (10) having minimum error were obtained by using the arbitrary values of the time-dependent variables x(t), y(t) and z(t) shown in Table 4.In this table, x(t) indicates the number of teachers exposed to mobbing, y(t) indicates the number of teachers who are not exposed to mobbing or gained resistance to it, and z(t) indicates the number of superiors who apply mobbing.A series of operations summarized respectively as following were carried out.
• The system proposed in Eq. ( 8) is expressed with ODEs.
• The ODE model was solved by RungeKutta45 with the Matlab R2023a, and then the values of parameters closest to the values in Table 4 (giving the minimum error) were found by using the lsqcurvefit function.• The values of parameters obtained by the least squares method were rewritten into the fractionalorder model suggested in Eq. ( 8).• The derivative order in Eq. ( 8), whose parameters have been already determined, was taken into account for α = 1, 0.99, 0.98, . . .and the derivative-order given the values closest to the values in Table 4 was calculated.The reason for reducing the derivative order is to investigate the delay in the system because there is a delay inherent in fractional derivatives.

Determining the model parameters
Annual study calendar in private education courses in Türkiye; it is planned to be at least eighteen and at most thirty-six weeks during the academic year [38].In this context, the variable values to be used for numerical calculations are given in Table 4.As stated earlier, the values x(t), y(t) and z(t) indicate the proportion of each compartment in the total population.If the model suggested in Eq. ( 8) is taken into account through ODEs, the system, η , is obtained.Equations in (23) are defined by positive initial conditions x(t 0 ) = x 0 , y(t 0 ) = y 0 and z(t 0 ) = z 0 for t 0 ≥ 0. The parameter values of Eqs. ( 23) have been investigated using the least squares method to obtain values close to the values in Table 4.For this, Eq. ( 23) has been solved with RungeKutta45 and suitable parameters are searched with the lsqcurvefit function.The Matlab codes used are located in Appendix A. By running the run_parameter_determining.m file, the warning "Local minimum possible.Its initial value is less than the default value of the function tolerance." is obtained.It has been lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the default value of the function tolerance.In this sense, its stopping criteria are the follows: optimization stopped because the relative sum of squares (r) is changing by less than options.Observe the following results: Relative change r = 1.78e − 08, TolFun = 1e − 06 (default).As a result, we have "Rate Constants" given as:

Selection of the derivative order that minimizes the error
In this section, the order of derivative, which minimizes the error, is searched according to Table 4.
The proposed fractional-order model in the Caputo sense in Eq. ( 8) with the parameter values showed in Eq. ( 24) are where x(0) = 0.741, y(0) = 0.021, z(0) = 0.004.The solutions have been found by selecting the derivative-order of (25) as 1, 0.99, and 0.98.Considering Table 4, the distances in R 3 between the variable values in this table and the values obtained from the solutions at the relevant time have been found.The distances for each row show how much they deviate from the true value.Then, to compare the derivative orders, the distances have been summed for each derivative order.
Since it has been observed that the sum of the distances after 0.98 gradually increased, the smaller derivative orders have not been given for the comparison.Relevant Matlab codes are provided in Appendix B. Related graphics can be seen in Figure 2, Figure 3, Figure 4, Figure 5.   4.
The sum of the calculated distances of the solution curves to the values in Table 4 for each derivative-order are as follows: Total_d_for_alpha_1 = 0.0566, Total_d_for_alpha_.99 = 0.1094, Total_d_for_alpha_.98 = 0.1947.
According to these results, the most appropriate derivative-order for the system in Eq. ( 25) is α= 1 as seen in Figure 2, Figure 3, Figure 4, Figure 5.In Figure 5, the relation between the real values and the solution curves of the system in Eq. ( 25) for each derivative-order are displayed graphically in R 3 .

Stable equilibrium point according to Table 4
Here, the calculations have been made for the equilibrium points of Eq. ( 8) according to the parameter values given in Eq. (24).In this context, the equilibrium points of system (8) are obtained as E 0 (3.0752, 0.44935, 0) and E 1 (0.31908, 0.209956, 0.03198).For these points, E 1 , the positive equilibrium point, is LAS.In other words, the solution of the system (8) converges to this point over time with initial conditions within a certain range.This situation can be clearly observed in Figure 6.25) for α = 1, 0.99, 0.98 and their comparison with values in Table 4.   25) for α = 1, 0.99, 0.98 and their comparison with values in Table 4.

Results and discussion
Mobbing, which can be defined as oppressive management and psychological violence carried out knowingly and in a certain way, such as disturbing, humiliating, excluding, ignoring, etc. for the purposes of pacifying, wearing out, and intimidating the employee, is a major problem in many public or private organizations.Studies have shown that people exposed to mobbing practices experience health problems such as depression, anxiety, stress disorder, high risk of suicide, insomnia, psychosomatic health complaints, musculoskeletal disorders, and fibromyalgia   (soft tissue rheumatism).Therefore, this problem can result in many problems that will negatively affect employee performance.Many individuals who encounter this problem receive psychological or medical support.Additionally, institutions take various precautions to minimize the possibility of their staff being subjected to mobbing.The fact that communities of individuals with various characteristics always interact cannot be ignored.In a given community, there is an interactive relationship between the compartments that are exposed to mobbing, those who are not exposed to mobbing, and those who practice mobbing (especially those at the management level).
The relationship between such compartments can be examined mathematically with differential equations based on the SIR model.In this way, time-dependent changes in compartment sizes can be predicted, and various control strategies can be presented as a result.
It is a fact that mathematical modeling is applied to many different fields today, from social sciences to medical sciences and engineering.The use of fractional order differential equations in mathematical modeling means that any delay in the system can be examined.In this study, the proposed mathematical model is presented with fractional-order differential equations, and the relationship among compartments based on mobbing is examined.The model has two different equilibrium points: the equilibrium point where individuals do not practice mobbing and the positive equilibrium point.From these points, the stability of the positive equilibrium point is shown.According to the analysis results, if there are individuals who practice mobbing in a closed working environment, mobbing will never disappear.With positive compartment numbers, the scenario continues permanently.In private or public institutions diagnosed with the presence of mobbing, our model provides the fact that intervention should be made without delay.
The order of derivatives and the parameters used in the model were analyzed in the graphical studies section.Here, the parameter values in sample A College were determined by using the time-dependent values of the variables representing the by x(t) the teachers who were exposed to mobbing, by y(t) the teachers who were not exposed to mobbing or gained resistance to mobbing and by z(t) the administrators who applied mobbing.Let us consider the evaluation to be sampled for N = 1000.These parameter values were written into the proposed model and the derivative-order was obtained as 1, which was closest to the data of A College.The error term was investigated with the derivative order being 1, 0.99, and 0.98.Here, the reason why the derivative order was examined by giving lower values is that it was desired to investigate any delay in the system.In this comparison, the concept of distance in R 3 was used by taking into account the equivalents of each point at each derivative order for A College.We have assumed that the total number of individuals is 1000 and the initial conditions of the proposed model are x(t 0 ) = 741 y(t 0 ) = 21 and z(t 0 ) = 4.These variables will be close to the values of x(t) = 319 y(t) = 210 and z(t) = 32 over time.While it was observed that there was a serious decrease in the number of teachers exposed to mobbing, it was concluded that there was a serious increase in the number of teachers who were not exposed to mobbing or who became resistant to mobbing and the number of administrators who practice mobbing.
Another important point in the model is the process of determining the parameters by using data to sample College A. Here, a Matlab 2023 program was used, and the relevant codes are included in Appendix A and Appendix B. Here, given the time-dependent changes in the number of individuals in any community with the proposed approach, this method Future predictions can be made easily using and so, a preliminary warning will be created for the development of control measures and planning strategies for mobbing.
In future studies, the parameters that would minimize the error between the real data will be determined by mathematical modeling with fractional-order differential equations containing different derivative orders in the system.The desire to be examined is the desire to examine whether there may be a delay in a particular variable, different from all variable values.

Figure 5 .
Figure 5.The solution curves for each variable of Eq. (25) for α = 1, 0.99, 0.98 and their comparison with values in Table4

Table 1 .
State variables in the proposed mathematical model.

Table 2 .
The parameters used in the proposed model and their meanings.

Table 3 .
Equilibrium points and LAS conditions.

Table 4 .
x(t), y(t) and z(t) values for college A.

Table 4 is
A_college_data.xls file in numerical calculations