Bullying is a pervasive issue that affects individuals across ages, genders, and cultures. It includes the intentional and persistent use of physical force, emotional or psychological aggressiveness to injure, threaten, or dominate another individual. While bullying has always existed, its prevalence and severity have increased in recent years due to the widespread use of technology and social media. This has become a common problem in schools (institutions of learning at different levels) and communities that has negative effects on the victims and society as a whole. Bullying can take different forms such as physical, verbal, relational, and cyberbullying. This literature review aims to examine the current state of research on bullying, its prevalence, effects on victims, and interventions to prevent and address bullying.

Bullying can lead the victim to commit crime, start engaging in drugs and smoking, and has psychological effects. These social vices have been studied extensively in the literature. In Nwajeri et al. [1], the authors studied the co-dynamic of drug trafficking and money laundering using fractional derivatives and also in Ngungu et al. [2], authors considered smoking in different age groups using fractal fractional derivatives and how government can intervene to curb the menace. In Alaka et al. [3], authors considered the psychological effect of the Coronavirus disease outbreak in a given population.

Bullying is a complex phenomenon that can be difficult to define, as it involves a wide range of behaviors and motivations. Some definitions emphasize the repetitive nature of bullying, while others focus on the power imbalance between the bully and the victim (or bullied). According to the Centers for Disease Control and Prevention (CDC), bullying is characterized as persistent, unwelcome aggressive behavior from one young individual or a group of youths who are not related or in current dating relationships. This behavior entails a perceived or observable power imbalance and occurs repeatedly or is highly anticipated [4].

For the purpose of the model under consideration, we divide the total student population into four compartments: namely the students with lack of family education X₁, that is, students who lack family education and not involving bullying; students with family education X₂, that is, student who have family education and not involving bullying; the bully students with lack of family education Y₁; the bully students with family education Y₂. The following nonlinear ordinary differential equations are generated by considering the interplay between the compartments:

where the total population at time is

The corresponding non-negative initial conditions to the Equation 1 are given below

The recruitment rates of school enrollment in the Equation 1 is denoted by A. m denote the rate of students who lack education from their families while ξ₁ and ξ₂ are progression to X₂ from X₁ and to X₁ from X₂, respectively. The rate of X₁ and X₂ graduate from schools are denoted as λ₁ and λ₂, respectively. The non-negative parameters μ₁ and μ₂ measure students who leave school due to the harassment of the bullies. The rates at which Y₁ and Y₂ bullies on other students are measured by β₁ and β₂, respectively. The rate of Y₁ and Y₂ who are expelled from school are denote by ν₁ and ν₂. Parameters ψ and ϕ measure the proportion of Y₂ who become Y₁ and proportion of Y₁ who become Y₂, respectively.

The key parameters associated with family education are:

The inclusion of family education is built on several underlying assumptions:

Tilahun et al. [45] discuss the power law correlation function and its implications for a time-dependent kernel. This discussion allows us to integrate the Caputo fractional order derivative into our model. As a result, the ordinary differential Equation 1 is reformulated into a Caputo fractional-order system as shown below:

where t∈[0, η], T∈ℝ, and Dtα0C denotes the Caputo fractional derivative of order α, where the memory index is denoted as α, 0 < α ≤ 1. Note that when α = 1 the fractional-order model reduces to the ordinary differential model (Equation 1). It is commonly known that Caputo's derivative offers greater dependability and adaptability in analytical applications. Due to initial condition properties, which are more physically interpretable for most problems, many researchers consider the Caputo operator over other derivatives. Consequently, we represented the dynamics of the bullying using fractional theory. Table 1 describe and defines the system parameters.

The section details various key definitions, lemmas, and ideas that are required to comprehend our proposed model.

Definition 2.1 [46, 47] Given an integrable function u:ℝ⁺ → ℝ, and α∈(n−1, n), n∈N. The left Caputo fractional derivative of order α of the function u is defined as

Definition 2.2 [47] The appropriate Riemann-Liouville fractional integral for the power-law kernel is given as

Lemma 2.3 [47] Assuming there is a function u(t)∈C([c, d], ℝ) of order α∈(0, 1). The existence-uniqueness of the fractional differential equation

is continuous and bounded and satisfies a Lipschitz condition.

This section will present solutions for the existence and uniqueness of the Caputo model for the Equation 2. Suppose that there exists a set of real valued Banach space functions denoted G on the interval B=[0,T] and let denote M=G(B)×G(B)×G(B)×G(B) with norm ||(X₁, X₂, Y₁, Y₂)|| = ||X₁||+||X₂||+||Y₁||+||Y₂||, where ||X1||=supt∈B|X1|,||X2||=supt∈B|X2|,||Y1||=supt∈B|Y1|,||Y2||=supt∈B|Y2|. We reformulate the right-hand side of Equation 2 as follows:

From understanding in Equation 3, our proposed Caputo fractional Equation 2 can be expressed as

where

and

Using the ideas in Lemma 2.3, Equation 4 yields

From here, our aim is to find the existence and uniqueness solutions, we assume the function vector

as:

(H1) ∃ ΩB,ΘB such that

(H2) ∃ LB>0 such that if D,D~∈M, then

Theorem 4.1. Under the continuity of B together with (H1), Equation 4 has at least one solution.

Proof . We shall arrive at the required conclusion using the Schauder fixed point theorem. Let us take a closed subset G of B as WB={D∈B:||D||≤ℝ,ℝ>0}, where WB is the operator defined as WB:G→G such that

which means that

From Equation 12, it follows that |Ω(D)|≤|D0|+TαΓ(α+1)[ΩB||D||+ΘB] and also WB∈G such that WB(D)⊂G. Also it reveals that the operator WB(D) is bounded. For completely continuity we proceed as follows. Let ϑ2<ϑ1∈[0,T] such that

Basically, we can see from Equation 13 that |WB(D)(ϑ2)|-|WB(D)(ϑ1)|→0 as ϑ₂ → ϑ₁. Hence, WB is an equicontinuous operator. With the help of the Arzela Ascoli theorem, we know that function WB is a completely continuous function and uniformly bounded. Again, by Schauder's fixed point theorem, we conclude that our proposed system has at least one solution.

Theorem 4.2. Suppose that (H2) holds. Then the considered Equation 1 has a unique solution if

Proof . If D,D~∈G, then

Hence, ||WB(D)-WB(D~)||≤TαΓ(α+1)LB||D-D~||, which completes the proof: from the contraction principle, the operator has a unique fixed point and, consequently, our proposed model has a unique solution.

In this section we examine the stability results of our proposed system. Hyers-Ulam (HU) stability is a significant concept in the analysis fractional equations, which includes models like the one for bullying dynamics among students. It pertains to the stability of functional equations when subjected to small perturbations, indicating that if a function nearly satisfies a functional equation, there exists a precise solution near this approximate one. This principle can be utilized for ordinary differential equations (ODEs) to guarantee that small alterations in initial conditions. By proving that the model is Hyers-Ulam stable, the study underlines the suitability of fractional-order modeling for social behavior systems like bullying, where the past influences present behavior. This is crucial for justifying the choice of fractional calculus in the proposed model, ensuring that the memory aspect contributes to a robust representation of the system dynamics. The stability analysis also implies that interventions based on the model's predictions can be implemented with the understanding that minor variations in their execution will not drastically affect the overall outcomes. This is important for policymakers and educators who need to rely on the model for decision-making, as it implies a degree of resilience to unexpected variations in program implementation or conditions. For stability of our proposed model we consider the following Lemma 3 and Theorem 4:

Lemma 3. The function D(ϑ) is a solution of

satisfies the given relation

Proof. The proof is very common in literature so we omit it, see Ngungu et al. [2].

Theorem 4. The Caputo fractional model under consideration is HU stable, if there exists

Proof. Let D∈G be any solution and D~∈G be a unique solution, then

Thus,

Hence, we conclude that the Caputo fractional model under consideration is HU stable.

Here, using Adams-Bashforth method, we provide the numerical algorithms for the Caputo fractional model under consideration. The formula for the Cauchy problem of the Caputo derivative is

From principle of the Caputo integral, Equation 19 is rewrite as

At the point ϑ_p+1 = (p+1)h and ϑ_p = ph, p = 0, 1, 2, 3, 4, ..., with h being the time step, Equation 20 will then be rewrite as:

which is the same as

For simplicity we take the right hand side of Equation 22 and using the Lagrange polynomial, Equation 22 can now be written as

where Equation 23 can be written as

For simplicity we evaluate the integrals in Equation 24, hence we have:

For simple substitution of Equations 25, 26 into Equation 24, we get the numerical algorithm for the Caputo derivative is as follows

Thus, in terms of our Caputo fractional social epidemiological model for the bullying, we get:

The diagram of our proposed model under the Caputo fractional operator is depicted in Figure 1, using the Adams-Bashforth method and taking into account the following initial conditions: X₁(0) = 500, X₂(0) = 110, Y₁(0) = 35, Y₂(0) = 20 with parameter values in Table 2. The goal of this section is to investigate the influence of changes in the values of these input parameters on system behavior. Our goal in studying many scenarios is to obtain a thorough grasp of how the system responds to changing conditions. By examining these numerical results, we hope to highlight the system's robustness and sensitivity to varied input variations. This section is critical in demonstrating the system's functioning under diverse situations, which deepens understanding of its responses and provides significant insights for scientific applications.

In the first simulation presented in Figures 1A, B, we noticed that the behavior of students who lack and with family education and not involving bullying under different values of fractional order. The trajectory changes with the fractional order, suggesting that the system's sensitivity to initial conditions or recent changes decreases as fractional order reduce from 1. Essentially, a lower fractional order makes the system less reactive to short-term fluctuations and more stable.

In Figures 1C, D, we illustrated the solution pathways of students with/without parental control and bullying other students of the system. In this simulation, we observed from Figure 1C that a fractional order is associated with a increase in bullying, unlike Figure 1D where a lower fractional order is associated with a decrease in bullying, these trajectories imply that interventions that have a more prolonged effect or that contribute to the “memory” of the system such as consistent policy enforcement or cultural shifts against bullying could be more effective than short-term solutions.

In Figures 2, 3, the effects of varying different parameters associated with family education are analyzed, providing insights into how changes in family dynamics and educational background influence the prevalence of bullying over time. Figures 2A, B present the dynamics of students without and with family education respectively, showing how the parameter variations (ξ₁ and ξ₂)impact students' stability. These parameters represent the rates of transition between having and lacking family education. The pink trajectory represents higher transition rates, indicating quicker movement from lacking family education to gaining it and vice versa. In both cases, the pink curves stabilize at higher values compared to the blue curve, indicating quicker movement from lacking family education to gaining it and vice versa. In both cases, the pink curves stabilize at higher values compared to the blue curves. This suggests that policies that facilitate family engagement or provide external support to simulate family education (e.g., community programs or school counseling) can lead to better stability and outcomes in students' behaviors. Improving family education rates will potentially reduce students' vulnerability to adverse behaviors by enhancing resilience against bullying. Figures 2C, D demonstrate the dynamics of bullying among students lacking family education and those with family education, respectively. The dynamics show that a higher rate of transition into a supportive environment leads to reduced numbers of bullying incidents over time. These results highlight the importance of family involvement as a preventative factor against bullying. Schools could collaborate with parents to strengthen family education programs, provide workshops on parenting skills, and offer resources aimed at creating supportive home environments. Furthermore, the transition rates (ξ₁ and ξ₂) could be influenced by school-led initiatives that directly involve parents in students' lives.

In Figures 3A, B depict the dynamics of students with and without family education for varying values of m which represents the rate of students lacking family education. When m increases from 0.05 to 0.09, it can be observed that the number of students without family education increases, as seen in the red curves. Conversely, the number of students with family education decreases, implying a growing population of students who lack sufficient support from their families. These figures show that increasing the rate of students without family education negatively impacts the student population's dynamics, increasing vulnerability. In real-world settings, this points to the need for government and school authorities to intervene with targeted programs for families that struggle to provide sufficient education and support. Policies aimed at subsidizing educational resources or counseling services could reduce m and improve student outcomes. Figures 3C, D show the dynamics of bullying incidents among students with and without family education as the rate m changes. With increased m the number of bullying incidents among students without family education also rises. This relationship indicates that when more students lack family education, the likelihood of involvement in bullying, either as a perpetrator or as a victim, also increases. For this results, schools and educational policymakers should focus on interventions that provide additional support to students who lack family education, such as mentorship programs, peer support initiatives, and after-school activities. These can serve as compensatory measures for students who are at higher risk due to a lack of family support. Providing targeted emotional and social learning initiatives in schools can also help mitigate bullying behaviors among students who may be predisposed to such behavior due to their family background.

In our ever more interconnected world, the capacity to model and anticipate the spread of bullying is essential for global readiness. Utilizing mathematical models, early identification and containment strategies can be guided, playing a crucial role in averting the swift global dissemination of bullying. Understanding the dynamics of a bullying model can aid in communicating the importance of preventive measures to the public. When people grasp how a bullying spreads among students, they are more likely to adhere to recommended behaviors like parental education and other preventive polices.

The students with lack of family education X₁, that is, students who lack family education and not involving bullying; students with family education X₂, that is, student who have family education and not involving bullying; the bully students with lack of family education Y₁; the bully students with family education Y₂.

Bullying has non-positive impacts on the mental-physical health, academic outcome, and social relationships of the victims. Those who are bullied face higher rates of sadness, anxiety, low self-esteem, and thoughts of suicide and post-traumatic stress disorders (PTSD) sometimes leading to mental health pathologies in the victims, prompting the need for a comprehensive examination of its cause and dynamics. This research has investigated the dynamics of bullying by introducing a parental education and has employed the fractional Caputo operator to model the system. The fractional Caputo model has been derived, and key findings related to this fractional model have been presented. Notably, stability analysis has been conducted, and the uniqueness and existence of the solution for the bully system have been examined. A novel numerical approach has been introduced to explore the system of bullying. The use of different orders of derivatives has revealed that the hypothesized model with the Caputo derivative has yielded more accurate and efficient graphical results. From the outcome of our result numerically, we have been able to show interventions have effect on reducing the prevalence of bullying; such interventions could campaign or policy shifts which gives better result than [44]. The outcomes of this study are expected to have provided valuable insights for policymakers, parents, and school authorities in their efforts to combat the spread of bullying in their communities or schools. Additionally, future iterations of this model may incorporate time-dependent control strategies based on optimal control theory to enhance bullying management and also test the model with real data by using the non-spatial form of the model developed at small spatial scales to understand the effect of bullying at different district schools in places like the USA at the census tract level.