Nonlocal Models in Biology and Life Sciences:
Sources, Developments, and Applications

Homogeneous and heterogeneous systems are fundamental concepts in mathematical biology that play crucial roles in understanding various biological phenomena. Homogeneous systems are the baseline for understanding the fundamental dynamics of biological systems before incorporating more complex factors. For instance, patterns are observed in nature and have visible regularities. These patterns recur in different contexts and can sometimes be modelled mathematically. According to Turing [3], “when two chemical species with different diffusion rates react with each other, the spatially homogeneous state may become unstable, thereby leading to a nontrivial spatial structure”. This looks counter-intuitive as it is expected to be a uniform distribution of the chemical substances through the diffusion term. In general, the diffusion term in RD equations contributes to transferring the concentrations throughout the spatial domain, and the linear or nonlinear reaction terms come from the mass action laws [4, 62]. Without loss of generality, we can consider a general form of RD equations with $m$-substances $\mathbf{u}(\mathbf{x},t)=(u_{1}(\mathbf{x},t),\ldots,u_{m}(\mathbf{x},t))^{T}$ at time $t>0$, and at $\mathbf{x}$ in a spatial domain $\Omega\subseteq\mathbb{R}^{n}$ as [4]:

$${}\frac{\partial\mathbf{u}}{\partial t}=\nabla\cdot(\mathbf{D}\nabla\mathbf{u})+\mathbf{f}(\mathbf{u}),$$

(1)

with non-negative initial conditions $\mathbf{u}(\mathbf{x},0)=\mathbf{u}_{0}(\mathbf{x})$ and appropriate boundary conditions, e.g., Dirichlet, no-flux, periodic, etc. Different types of boundary conditions have been used in RD equations, depending on the nature of the problems; no-flux and periodic boundary conditions are among the most familiar. Here, $\mathbf{D}$ is a diagonal matrix called the matrix of diffusivity. The term $\nabla\cdot(\mathbf{D}\nabla\mathbf{u})$ is the self-diffusion of the substances, taking care of the rate of change of the concentrations at a particular point in the spatial domain. The matrix $\mathbf{D}$ is not a diagonal for the case of cross-diffusion. The initial conditions for all the species are assumed to be non-negative since the concentration is always non-negative. The vector-valued function $\mathbf{f}(\mathbf{u})$ represents the reactions between the species in the spatial domain $\Omega$. These reactions could be the variables’ source, sinks, and interactions. Therefore, modelling heterogeneous systems is essential for understanding the impact of individual variability on population dynamics, disease spread, and other biological processes.

Depending on the parameters and interactions between organisms, the RD equation (1) exhibits different types of dynamic non-homogeneous solutions, such as travelling waves, wave-trains, spatio-temporal chaos, etc. In addition, reaction-diffusion systems can produce non-homogeneous stationary patterns such as spots, stripes, and target patterns [4]. In 1937, Fisher and Kolmogorov, Petrovskii, and Piskunov introduced the scalar RD equation [63]:

$${}\frac{\partial u}{\partial t}=d\frac{\partial^{2}u}{\partial x^{2}}+g(u),$$

(2)

and studied the existence of the travelling wave solutions along with the propagation of the wave speed. The equation (2) is known as the Fisher-KPP equation. Researchers have been studying and applying this Fisher-KPP model in various fields, e.g., physics, chemistry, biology, and medicine. A particular nonlinear form $g(u)=u(1-u)$ of the RD equation (2) has been taken a lot of attention in capturing different phenomena through mathematical modelling, e.g., spatial spreading of cell populations [64], wound-healing cell migration [65, 66], protein misfolding [67, 68], biological invasion [69, 70, 71], etc.

As discussed earlier, the RD equations (1) are capable of producing different types of stationary and non-stationary non-homogeneous patterns. Several mechanisms are available for getting non-stationary patterns [4]. Turing suggested a mechanism to obtain some of the non-homogeneous stationary patterns mathematically, called the Turing pattern [3]. One of the main conditions for existing Turing patterns requires at least two interacting populations for the RD equations. Therefore, the usual Turing mechanism can not be applied to the single-species local RD model (2). However, it can be applied to a single species RD model with nonlocal interactions. We first describe a nonlocal version of the Fisher-KPP RD model for a single species population (2), then move to the multi-species populations and other nonlocal interactions.

Without loss of our generality, we assume that $u$ is a dimensionless population density, and it is defined in the one-dimensional spatial domain $\Omega=[-L,L]\subset\mathbb{R}$, $L>0$. In the considered case $g(u)=u(1-u)$, the population has the reproduction term, which is proportional to its density and available resources $(1-u)$ [62, 72]. For the local interaction, the population’s consumption rate at any given spatial point $x\in\Omega$ is proportional to its population’s density at the same spatial point with the total amount of available resources unity. In general, the individual located at a spatial point $x$ accesses the resources at its location and some in its neighbourhood locations [41, 72, 73, 74, 75, 76]. We assume that the individuals located at the spatial point $y$ access the resources of the point $x$ with the probability density function $\phi(x-y)$. In this case, the RD equation (2) modified into an integro-differential equation

$${}\frac{\partial u}{\partial t}=d\frac{\partial^{2}u}{\partial x^{2}}+u(1-\phi\ast u),$$

(3)

with non-negative initial conditions and periodic boundary conditions over the domain $[-L,L]$. Here $\phi\ast u$ is the convolution of $u$, and we use $`\ast$’ as the convolution operator for the rest of this manuscript. The convolution $\phi\ast u$ in one spatial dimension is defined by

$$(\phi\ast u)(x,t)=\int_{-\infty}^{\infty}\phi(x-y)u(y,t)dy.$$

In the RD equation (2), the diffusion term contributes to the net density change in the spatial locations and is much slower than the involvement of the consumption resources. But in the nonlocal model (3), the convolution term captures the long-range interaction in the spatial domain, one of the key factors missing in the local model. The kernel function $\phi$ satisfies the normality condition to share the same homogeneous steady-state solution corresponding to the local model, and the condition is

$${}\int_{-\infty}^{\infty}\phi(y)dy=1.$$

The nonlocal model (3) has two homogeneous stationary solutions $u=0$ and $u=1$. In the case of the local model, the solution $u=1$ is stable under heterogeneous perturbations, but it may not be stable for the nonlocal model [62]. A brief analysis of this single species nonlocal model is given in the next section.

Researchers have studied nonlocal RD models for the triangular, parabolic kernel functions to account for different biological phenomena [77, 78, 79, 80, 46, 81, 82]. These kernel functions are known as the truncated kernel. Such non-truncated kernel functions have been the subject of intensive studies. For instance, Britton introduced an advanced version of the single species model (3) as:

$${}\frac{\partial u}{\partial t}=d\Delta u+u(1+\alpha u-(1+\alpha)\psi\ast u),$$

(4)

where $\alpha u$ ($\alpha>0$) represents the local aggregation and $-(1+\alpha)\psi\ast u$ is the disadvantage in the global population as the resources may go into depletion of being high population density [83]. The author studied the nonlocal model with the Laplacian kernel function [46]. Sometimes, species take time to move; in this case, the correct average population density is a spatio-temporal weighted average toward the current time and position. Taking this spatio-temporal factor into account, Britton introduced a space-time nonlocality in the single species model (4) as [84]:

$${}\frac{\partial u}{\partial t}=d\Delta u+u(1+\alpha u-(1+\alpha)\psi\ast\ast\leavevmode\nobreak\ u),$$

(5)

where the double asterisk represents the convolution in space and time and is given by

$$(\psi\ast\ast\leavevmode\nobreak\ u)(\mathbf{x},t)=\int_{-\infty}^{t}\int_{\Omega}\psi(\mathbf{x}-\mathbf{y},t-s)u(\mathbf{y},t)d\mathbf{y}ds.$$

Likewise, the spatial kernel function $\Phi$ in (3), the spatio-temporal kernel function $\psi$ satisfies some conditions mentioned in [84, 85]. In addition, these types of nonlocal models with non-truncated kernels can be converted into coupled RD equations [83, 84]. A brief formulation of this conversion and its applications are given in the upcoming section.

In cell biology, nonlocal interactions are interactions between cellular components that extend beyond immediate neighbours and may involve long-distance signalling or communication. These interactions are critical in cellular processes and behaviours such as cell migration, pattern formation, and collective cell responses. For example, tumors are abnormal masses of tissue consisting of cells with different phenotypes (e.g., cancer stem cells) [86, 87, 88], and they are heterogeneously distributed in the body. Motivated by the paradoxical findings of an agent-based model [89], Hillen et al. have proposed a mathematical model [90]:

	$\displaystyle\frac{\partial u(x,t)}{\partial t}$	$\displaystyle=D_{u}\Delta u+\delta\gamma\int_{\Omega}\psi(x,y,p(x,t))u(y,t)dy,$		(6)
	$\displaystyle\frac{\partial v(x,t)}{\partial t}$	$\displaystyle=D_{v}\Delta v+(1-\delta)\gamma\int_{\Omega}\psi(x,y,p(x,t))u(y,t)dy-\alpha v+\rho\int_{\Omega}\psi(x,y,p(x,t))v(y,t)dy,$

where $u(x,t)$ and $v(x,t)$ are the densities of cancer stem cells and non-stem cancer cells at the location $x$ and time $t$. Here, $p(x,t)=u(x,t)+v(x,t)$ is the total tumour density, and $\gamma$ and $\rho$ are the number of cell cycles times per unit time of cancer stem cells and non-stem cancer cells, respectively. The parameter $\delta$ describes the fraction of cancer stem cell divisions, and $\alpha$ is the tumour cell death rate. $D_{u}$ and $D_{v}$ are the usual diffusion coefficients of cancer stem cells and non-stem cancer cells, respectively. The kernel function $\psi(x,y,p)$ describes the rate of progeny contribution to the location $x$ from a cell at location $y$, per “cell cycle time” and satisfies the normality condition

$$\int_{\Omega}\psi(x,y,p(x,t))dx\leq 1.$$

Gene silencing and RNA interference are exciting fields of study that have promising implications in health and industry [91]. For instance, a synthetic RNA molecule is designed to selectively target and mute the production of one particular gene in small interfering RNA technology. When it enters cells, it attaches itself to the target mRNA, causing it to be degraded and stopping the creation of new proteins. There are some more key technologies and models associated with gene silencing and RNA interference: short hairpin RNA [92], micro RNA [93], CRISPR/Cas9 technology [94], Antisense Oligonucleotides [95], Dicer Knockout models [96], etc. It is essential to model genetic switches and genetic networks in order to comprehend the complex regulatory processes that exist within cells. The dynamics of genetic switches and networks have been simulated and examined using a variety of computational and mathematical models. Some common approaches are Boolean networks, ordinary differential equations, stochastic models, agent-based models, and machine learning approaches. Furthermore, nonlocal models are important in various areas beyond DNA, RNA, and protein modelling. The use of microarray technology for placing the full human genome on chips is one such application.

A multi-scale modelling technique is required to represent a cell as a composite structure incorporating microtubules and organelles. Several models have been created to simulate these components’ interactions and dynamics, such as continuum mechanics models [97, 98, 99]. Understanding the methods by which a cell develops, duplicates its components, and divides into two daughter cells is fundamental to the study of cell cycles. Modelling cell cycles with fractional derivatives presents a mathematical technique that reflects the nonlocal and memory-dependent character of cellular dynamics. In addition, rheology and hysteresis are more frequently linked with materials and physical systems, and they may also be used to describe some aspects of cellular activity, including cell mechanics and cycles [100]. Understanding the interaction of grain shape, material characteristics, and nonlocal rheology is critical for forecasting and regulating the behaviour of dense granular flows [101].

Different types of biological rhythms occur in various physiological and behavioural processes in living organisms, and they are essential for maintaining homeostasis and optimizing physiological functions in organisms. In addition, angiogenesis is the physiological process that results in the formation of new blood vessels from pre-existing vessels [102]. This complex mechanism is required for a variety of physiological activities, including embryonic development and wound healing. However, angiogenesis has been linked to a number of pathological disorders, including cancer and inflammatory illnesses. Furthermore, biological cells are complex structures with distinct rheological properties, and understanding their mechanical responses is critical for a variety of applications, including disease diagnostics, medication delivery, and tissue creation. The study of these rhythms provides insights into how living organisms adapt to their environments and regulate various processes over time.

Nonlocal interactions on networks in biology refer to interactions between nodes in a network (such as genes, proteins, or cells) that extend beyond their direct neighbours or immediate connections. These interactions serve as essential for understanding the processes and characteristics of complex biological systems in a network. Furthermore, these interactions on networks are relevant in various biological systems, including signalling networks, ecological food webs, and neuronal networks.

Neurons in the brain communicate and interact through synapses, which can involve long-range connections. Nonlocal models have been employed to describe neurons’ collective activity and synchronization in neural networks. In addition, Alzheimer’s disease (AD) is also growing fast in the current days, and researchers have been trying to understand the mechanisms behind its propagation. Many factors are involved in AD progression, but two protein families, amyloid-beta and tau proteins, are known to be the main contributors [103, 104]. In [105], authors have considered the heterodimer model to incorporate the interaction between these two proteins, each consisting of healthy and toxic densities. Also, they derived the network model corresponding to the PDE model for applying it to AD progression in the brain connectome. Authors in [106] modified the heterodimer model by introducing nonlocal interactions into the model. They studied the following nonlocal network model

	$\displaystyle\frac{du_{j}}{dt}$	$\displaystyle=-\sum_{k=1}^{N}L_{jk}u_{k}+u_{j}(a_{0}-a_{1}u_{j})-\frac{a_{2}u_{j}}{1+c_{u}u_{j}}\Phi_{j}\ast\widetilde{u_{j}},$		(7a)
	$\displaystyle\frac{d\widetilde{u}_{j}}{dt}$	$\displaystyle=-\sum_{k=1}^{N}L_{jk}\widetilde{u}_{k}-\widetilde{a}_{1}\widetilde{u}_{j}+a_{2}\widetilde{u}_{j}\Phi_{j}\ast\bigg{(}\frac{u_{j}}{1+c_{u}u_{j}}\bigg{)},$		(7b)
	$\displaystyle\frac{dv_{j}}{dt}$	$\displaystyle=-\sum_{k=1}^{N}L_{jk}v_{k}+v_{j}(b_{0}-b_{1}v_{j})-b_{3}\widetilde{u}_{j}v_{j}\widetilde{v}_{j}-\frac{b_{2}v_{j}}{1+c_{v}v_{j}}\Phi_{j}\ast\widetilde{v}_{j},$		(7c)
	$\displaystyle\frac{d\widetilde{v}_{j}}{dt}$	$\displaystyle=-\sum_{k=1}^{N}L_{jk}\widetilde{v}_{k}-\widetilde{b}_{1}\widetilde{v}_{j}+b_{3}\widetilde{u}_{j}v_{j}\widetilde{v}_{j}+b_{2}\widetilde{v}_{j}\Phi_{j}\ast\bigg{(}\frac{v_{j}}{1+c_{v}v_{j}}\bigg{)},$		(7d)

where $u_{j}$ and $\widetilde{u}_{j}$ are the densities of healthy and toxic amyloid beta, and $v_{j}$ and $\widetilde{v}_{j}$ are the densities of healthy and toxic tau proteins at the node $j$ in a given graph. Here, $N$ is the total number of nodes and $L_{ij}$ are the entries of the Laplacian matrix for the given network. The convolution and the Laplacian of a graph are defined in [105, 106].

There are nonlocal mathematical models on neural networks, especially convolutional neural networks. The traditional neural network blocks feature representations in a local sense; however, long-range dependencies have significant practical learning problems. Here, we describe a nonlocal network in the context of traditional image classification tasks and make a comparison among different neural networks. Suppose $X=[X_{1},\ldots,X_{M}]$ is an input sample with $X_{j}(j=1,\dots,M)$ as the feature at position $j$. A nonlocal output signal at position $j$ can be defined as [107]:

$$Z_{j}=X_{j}+\frac{W_{Z}}{\mathcal{C}_{j}(X)}\sum_{\forall k}\omega(X_{j},X_{k})g(X_{k}),$$

where $W_{Z}$ is the weight matrix, $C_{j}(X)$ is the normalization factor, $\omega$ is a pairwise function which computes a scalar between the positions $j$ and $k$, and $g$ is a given function, e.g., $g(X_{k})=W_{g}X_{k}$. Incorporating these nonlocal blocks into a residual network, the network model can be written as:

$$Z^{n+1}:=Z^{n}+\mathcal{F}(Z^{n};W^{n}),\leavevmode\nobreak\ \leavevmode\nobreak\ n=0,1,\ldots,N,$$

where $W^{n}$ is a parameter set with $N$ number of network blocks with the initial input sample $Z^{0}=X$.

Genetic regulatory networks (GRNs) serve an important role in coordinating gene expression and orchestrating complex developmental processes in cells. These networks involve transcription factors binding to particular DNA sequences, cis-regulatory elements, and other substances, resulting in a dynamic and linked web of regulatory interactions. Understanding the abnormal alterations in GRNs gives insight into several disease pathways, such as cancer and developmental disorders [108]. Stochastic transcriptional regulation adds a layer of complexity to GRNs, altering cellular heterogeneity, phenotypic variability, and the general adaptability of biological systems [109, 110, 111]. Furthermore, pattern memory in GRNs is important in cellular decision-making processes, allowing cells to respond to the same stimuli differently based on their past experience.

The spatio-temporal dynamics of biomolecular networks and biochemical systems relate to the complex, dynamic patterns of molecular interactions that occur both in space (inside cellular compartments) and over time [112]. The spatial organization impacts many biomolecular interactions that occur in various cellular compartments, such as the nucleus, cytoplasm, and organelles. Furthermore, spatiotemporal dynamics in metabolic networks are critical for coordinating energy generation, nutrient use, and maintaining cellular homeostasis [113]. Nonetheless, spatial gradients of signalling molecules contribute to the generation of concentration gradients, impacting cell fate and pattern formation throughout development. The dynamics of transcriptional regulation involve both local and nonlocal components, reflecting the complex interplay of molecular processes within a cell as well as interactions that extend beyond immediate neighbours. Local dynamics often relate to nearby activities, such as interactions between genes and transcription factors. In contrast, nonlocal dynamics involve events that affect distant genomic areas or operate across distant locations and temporal dimensions.

Nonlocal dispersal in biology refers to the movement or migration of individuals, cells, or organisms over extended distances, often beyond their immediate vicinity. This phenomenon plays a crucial role in shaping populations’ distribution, dynamics, and persistence in various ecological and biological systems. Species expand their range of interactions through dispersal, which has many forms. Till now, we have considered the simplest type of dispersal as random diffusion, which captures the movement of organisms between adjacent spatial locations. A different type of dispersal that gives rise to nonlocality is studied in the literature, where organisms can travel some distance [114]. Periodic travelling wave solutions have been studied in prey-predator models with nonlocal dispersal to capture cyclic population dynamics. In [46], the author considered the Rosenzweig-MacArthur model with nonlocal dispersal as:

	$\displaystyle\frac{\partial u}{\partial t}$	$\displaystyle=\mathcal{P}u+u(1-u)-\frac{\alpha uv}{1+\alpha u},$		(8)
	$\displaystyle\frac{\partial v}{\partial t}$	$\displaystyle=\mathcal{P}v+\frac{\alpha uv}{\beta(1+\alpha u)}-\frac{v}{\beta\gamma},$

where

$$\mathcal{P}u=\frac{1}{l}\int_{-\infty}^{\infty}\phi\bigg{(}\frac{x-y}{l}\bigg{)}u(y)dy-u(x)$$

denotes the nonlocal dispersal corresponding to $u$ at a spatial point $x$ with $l$ as the spatial scale representing the dispersal distance [114]; $u(x,t)$ and $v(x,t)$ are the population densities of the prey and predator, respectively, at the spatial point $x$ and at time $t>0$. This type of nonlocal model is the derivative-free integral equation for diffusion [115]. To study the oscillatory dynamics of the model, the author converted this model into normal form near the Hopf bifurcation using the complex Ginzburg-Landau equation [116]. He estimated the wavelength and amplitude of the periodic travelling wave generated by the invasion of the coexisting steady-state for the parameter values close to the Hopf bifurcation threshold. Assuming specific properties for kernel function $\phi$ (non-negativity, being even, and normality), the nonlocal dispersal operator can be approximated for small $l$ as [114]:

$$\mathcal{P}u=\frac{l^{2}}{2}\bigg{(}\int_{-\infty}^{\infty}\phi(z)z^{2}dz\bigg{)}\frac{\partial^{2}u}{\partial x^{2}}+O(l^{3}).$$

On the other hand, for $l\rightarrow\infty$ and for a periodic $u$ with periodicity $L$, the operator $\mathcal{P}$ can be expressed as [114]:

$$\mathcal{P}u=\frac{1}{L}\int_{0}^{L}u(y)dy-u(x).$$

Fractional nonlocal models have gained significant attention in biology due to their ability to capture the memory effect, anomalous transport, and long-range interactions, which are common in biological systems. Fractional calculus has been employed to study drug release from polymeric matrices used in drug delivery systems. These models account for drug diffusion’s nonlocal and memory effects within the polymer, providing a more accurate representation of drug release kinetics. Different nonlocal RD models have been used in biology and other fields. Researchers have used fractional nonlocal models to understand cancer cell populations’ dynamics in cancer diseases. Before going to the fractional model, we first consider a modified nonlocal RD model (3) for nonlocal consumption of resources as [117]:

$${}\frac{\partial u}{\partial t}=d\Delta u+au^{m}(1-b\phi\ast u^{n})-cu.$$

(9)

The parameters $m$ and $n$ are positive integers, $m=1,n=1$ and $m=2,n=1$ correspond to the asexual and sexual reproductions, respectively. This model is an advanced version of the nonlocal model kinds of (3). Here, $a$, $b$ and $c$ are the parameters, and $\phi\ast u^{n}$ is the convolution. The last term in (9) represents the mortality with a rate constant $c$. In cancer disease, abnormal changes in cell density cause significant morbidity and mortality. Numerous mathematical models in the literature account for this disease propagation [118]. The fractional RD model is one of them, and many researchers have been focused on it. The time-fractional RD equation is used in the mathematical model to incorporate such an evolution process, and it is given by [117]:

$${}D^{\alpha}u=d\Delta u+au^{2}(1-b\phi\ast u)-cu,$$

(10)

where $D^{\alpha}u$ is the Caputo fractional derivative, and it is defined by

$$D^{\alpha}u(\mathbf{x},t)=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}\frac{u^{\prime}(\mathbf{x},s)}{(t-s)^{\alpha}}ds,\leavevmode\nobreak\ \leavevmode\nobreak\ 0<\alpha<1.$$

Here, $u^{\prime}(\mathbf{x},s)$ is the first-order partial derivative of $u(\mathbf{x},t)$ with respect to $t$ and evaluated at $t=s$. These models have many applications, and some of them are mentioned in the next section.

Nonlocal epidemiological models have been used to study the spread of infectious diseases in biological systems, considering interactions beyond immediate neighbours and capturing the effects of long-range dispersal and spatial dependencies. In 1927, based on certain assumptions, Kermack and McKendrick proposed a deterministic SIR model to capture a transmitted viral or bacterial agent in a closed population [14, 119]. Here, we first describe the SIR model and its generalization, incorporating the space into the model. Let $S(t)$ be the density of a population susceptible to a disease but not yet infected at the time $t\geq 0$. Suppose $A(\theta)$ is the expected infectivity of an individual who has been infected $\theta$ time ago. Then, the Kermack-McKendrick model follows the integral differential equation [14]:

$${}\frac{dS}{dt}=S(t)\int_{0}^{\infty}A(\theta)\frac{dS}{dt}(t-\theta)d\theta.$$

(11)

Based on the linearization of the integral differential equation (11), Kermack and McKendrick derived invasion criteria for which the disease invades individuals. To derive those criteria, we suppose that the whole population is susceptible at the epidemic’s beginning and $S(0)=S_{0}$. Also, we assume that the linearized solution at $S_{0}$ is of the form $c_{0}\exp(rt)$, and in this case, the characteristic equation is given by

$$1=S_{0}\int_{0}^{\infty}A(\theta)e^{r\theta}d\theta,$$

and the number

$$R_{0}=S_{0}\int_{0}^{\infty}A(\theta)d\theta$$

describes the epidemic’s growth in the initial phase of the disease propagation. The invasion criteria of the epidemic is $R_{0}$, which is greater than one. For the special choice $A(\theta)=\beta\exp(-\gamma\theta)$, where $\beta,\gamma>0$, the quantity

$$I(t)=-\frac{1}{\beta}\int_{0}^{\infty}A(\theta)\frac{dS}{dt}(t-\theta)d\theta$$

represents the number of infected individuals at the time $t$. If $R(t)$ is the number of individuals who recovered from the infection and may be infected again, then we can simplify all these equations in the following ODE system

	$\displaystyle\frac{dS}{dt}$	$\displaystyle=-\beta S(t)I(t),$		(12)
	$\displaystyle\frac{dI}{dt}$	$\displaystyle=\beta S(t)I(t)-\gamma I(t),$
	$\displaystyle\frac{dR}{dt}$	$\displaystyle=\gamma I(t).$

In 1965, Kendall extended this SIR model (12) into a nonlocal epidemic model [48]. Let $S(x,t)$, $I(x,t)$, and $R(x,t)$ be the densities of the susceptible, infected and removed individuals at the spatial location $x\in\mathbb{R}$ and time $t\geq 0$ such that their sum $S+I+R$ remains time-independent. We assume that the infected individuals are infectious and that the rate of infection is

$$\beta\int_{-\infty}^{\infty}I(y,t)\Psi(x-y)dy,$$

where $\beta$ is a positive constant and $\Psi$ is the normalized kernel function in $\mathbb{R}$. Taking these into the SIR model, the modified model becomes

	$\displaystyle\frac{\partial S}{\partial t}$	$\displaystyle=-\beta S(x,t)\int_{-\infty}^{\infty}I(y,t)\Psi(x-y)dy,$		(13)
	$\displaystyle\frac{\partial I}{\partial t}$	$\displaystyle=\beta S(x,t)\int_{-\infty}^{\infty}I(y,t)\Psi(x-y)dy-\gamma I(x,t),$
	$\displaystyle\frac{\partial R}{\partial t}$	$\displaystyle=\gamma I(x,t).$

Nonlocal epidemiological models have been applied in biology to study the spread of infectious diseases in spatially extended and heterogeneous environments, such as urban areas with varying population densities or regions with different levels of healthcare access [120]. These models incorporate long-range interactions and spatial dependencies, allowing for a more realistic representation of disease transmission dynamics.

The COVID-19 pandemic was a global health crisis caused by the SARS-CoV-2 virus in a heterogeneous environment [121]. This disease spreads rapidly through respiratory droplets when an infected person talks, coughs, or sneezes, and new variants of the virus continuously emerge, which causes concern among health experts around the world. Different types of measurements have been implemented to control the transmission of the virus, including targeted lockdowns, ramped-up testing, and widespread vaccination efforts. In terms of mathematical modelling, numerous mathematical frameworks have been proposed to investigate the transmission dynamics of the COVID-19 pandemic [122, 123, 124], and we highlight one of them where the spatial domain is considered. The novel coronavirus infects contacts during cold-chain transmission, so the study of the characteristics of COVID-19 dynamics can be more accurate in consideration of long-range population diffusion. Suppose $S(x,t)$, $V(x,t)$, $E(x,t)$, $I_{a}(x,t)$, and $I_{s}(x,t)$ are respectively for the populations of susceptible, fully vaccinated, latent infection, asymptomatic infected, and symptomatic infected at a point $x\in\mathbb{R}$ and time $t>0$, with $B(x,t)$ as the SARS-CoV-2 virus in the environment. Following [123], we consider the spatio-temporal nonlocal mathematical model for COVID-19 as:

	$\displaystyle\frac{\partial S(x,t)}{\partial t}$	$\displaystyle=d_{s}\mathcal{D}[S(x,t)]-(\mu(x)+\eta(x))S(x,t)+\Lambda(x)-S(x,t)f_{s}(x,t),$		(14)
	$\displaystyle\frac{\partial V(x,t)}{\partial t}$	$\displaystyle=d_{v}\mathcal{D}[V(x,t)]+\eta(x)S(x,t)-(\mu(x)+p(x))V(x,t)-V(x,t)f_{v}(x,t),$
	$\displaystyle\frac{\partial E(x,t)}{\partial t}$	$\displaystyle=d_{e}\mathcal{D}[E(x,t)]+V(x,t)f_{v}(x,t)+S(x,t)f_{s}(x,t)-(\mu(x)+\gamma(x))E(x,t),$
	$\displaystyle\frac{\partial I_{a}(x,t)}{\partial t}$	$\displaystyle=d_{a}\mathcal{D}[I_{a}(x,t)]+\epsilon(x)\gamma(x)E(x,t)-(\mu(x)+r_{a}(x)+\zeta_{a}(x))I_{a}(x,t),$
	$\displaystyle\frac{\partial I_{s}(x,t)}{\partial t}$	$\displaystyle=d_{s}\mathcal{D}[I_{s}(x,t)]+(1-\epsilon(x))\gamma(x)E(x,t)-(\mu(x)+r_{s}(x)+\zeta_{s}(x))I_{s}(x,t),$
	$\displaystyle\frac{\partial B(x,t)}{\partial t}$	$\displaystyle=d_{b}\mathcal{D}[B(x,t)]+\delta_{e}(x)E(x,t)+\delta_{a}(x)I_{a}(x,t)+\delta_{s}(x)I_{s}(x,t)-c_{s}(x)B(x,t),$

with non-negative initial conditions and Neumann boundary conditions. Here, $\mathcal{D}$ is the dispersal convolution operator defined by

$$\mathcal{D}[u(x,t)]=\int_{\Omega}\psi(x-y)[u(y)-u(x)]dy,$$

where $\psi(x-y)$ is the probability that an individual at position $y$ moves to position $x$ in the space $\Omega$. The functions $f_{s}(x,t)$ and $f_{v}(x,t)$ are defined as

$$f_{s}(x,t)=\int_{\Omega}[\beta_{e}(x,y)E(x,t)+\beta_{a}(x,y)I_{a}(x,t)+\beta_{s}(x,y)I_{s}(x,t)+\beta_{b}(x,y)B(x,t)]dy,$$

$$f_{v}(x,t)=\int_{\Omega}[\alpha_{e}(x,y)E(x,t)+\alpha_{a}(x,y)I_{a}(x,t)+\alpha_{s}(x,y)I_{s}(x,t)+\alpha_{b}(x,y)B(x,t)]dy,$$

where $\beta_{k}(k=e,a,s)$ are the rate of infection of infected people to susceptible people and $\alpha_{k}(k=e,a,s)$ are the rate of infection of infected people to those vaccinated, and $\beta_{b}$ and $\alpha_{b}$ are the rate of infection of the virus in the population. The other parameters used in the model (14) are given in [123]. The authors have shown that the effective rate of vaccination is one of the best measures to prevent and control the transmission of COVID-19 under current medical conditions.

Nonlocality in immunology refers to the phenomena in which immune system interactions and signalling processes extend beyond the immediate area of the immune cells involved. T-cell receptor (TCR) activation, for example, refers to the assumption that T-cell responses are not limited to the local contact between the TCR and its antigen. T-cell activation, on the other hand, is a networked and dynamic process that includes interactions at remote places [125]. It can lead to systemic effects, influencing immune responses and inflammation throughout the body, which has implications for autoimmune diseases and chronic inflammatory conditions [126]. Furthermore, nonlocal TCR activation is important in cancer immunotherapy because activated T cells may target cancer cells in distant metastatic locations. [127]. In T memory cell homeostasis, for example, the bone marrow serves as a niche and interacts with stromal cells, and cytokines contribute to their survival and nonlocal homeostasis [128]. These T cells grow and contract in distinct stages, including interactions and signals that reach beyond the T cells’ local surroundings. During T cell activation and proliferation, cytokines such as interleukins are generated not just locally at the site of activation but also systemically. These cytokines assist in the nonlocal regulation of T-cell proliferation [129].

Relaxation-time models are pretty interesting, especially in the context of life sciences. These models help us understand the dynamics of biological systems and provide insights into the time-dependent behaviour of various processes. For instance, the nonlocal model (5) in space and time can be considered as the relaxation-time nonlocal model. A similar type of extension can be applied to the model (13). We first include the spatial variable into the model (13) to account for the random movement of individuals as [130]:

	$\displaystyle\frac{\partial S}{\partial t}$	$\displaystyle=d_{1}\Delta S-\beta S(x,t)\int_{-\infty}^{\infty}I(y,t)\Psi(x-y)dy,$		(15)
	$\displaystyle\frac{\partial I}{\partial t}$	$\displaystyle=d_{2}\Delta I+\beta S(x,t)\int_{-\infty}^{\infty}I(y,t)\Psi(x-y)dy-\gamma I(x,t),$
	$\displaystyle\frac{\partial R}{\partial t}$	$\displaystyle=d_{3}\Delta R+\gamma I(x,t),$

subject to the Neumann boundary condition. Here $d_{1}$, $d_{2}$, and $d_{3}$ are the diffusion rates for the susceptible, infective and removed individuals. Studying a time delay to account for the latent period on the spread of the disease, we have to consider the following model:

	$\displaystyle\frac{\partial S}{\partial t}$	$\displaystyle=d_{1}\Delta S-\beta S(x,t)\int_{-\infty}^{t}\int_{-\infty}^{\infty}I(y,s)\Psi(x-y,t-s)dyds,$		(16)
	$\displaystyle\frac{\partial I}{\partial t}$	$\displaystyle=d_{2}\Delta I+\beta S(x,t)\int_{-\infty}^{t}\int_{-\infty}^{\infty}I(y,s)\Psi(x-y,t-s)dyds-\gamma I(x,t),$
	$\displaystyle\frac{\partial R}{\partial t}$	$\displaystyle=d_{3}\Delta R+\gamma I(x,t),$

where the kernel function describes the interaction between the infective and the susceptible individuals at location $x$ and the present $t$, which occurred at location $y$ and at earlier instance $s$. This kernel function satisfies some conditions which are mentioned in [119].

The nonlocal effects and time delay may come in the model when the mature death and diffusion rates are age-independent. To derive such a delay model, we follow the work [131] where the author considered a single-species population model with two age classes and a fixed maturation period living in a spatial transport field. Suppose $u(x,a,t)$ is the population density at a spatial point $x\in\Omega\subset\mathbb{R}$ and time $t\geq 0$ with age $a\geq 0$. Following [131], the single-species model given by

$${}\frac{\partial u}{\partial t}+\frac{\partial u}{\partial a}=\frac{\partial}{\partial x}\bigg{(}D(a)\frac{\partial u}{\partial x}+Bu\bigg{)}+d(a)u,$$

(17)

with the natural boundary condition. Here, $D(a)$ and $d(a)$ are the diffusion and death rates at the age $a$, respectively, and $B$ is the transport velocity of the field. For simplicity, we can assume that the population has two age stages: immature and mature. We also assume that $\tau\geq 0$ is the maturation time and $r$ is the life limit of an individual species. Then, the total matured population density is given by $\int_{\tau}^{r}u(x,a,t)da,\leavevmode\nobreak\ \leavevmode\nobreak\ x\in\Omega,\leavevmode\nobreak\ \leavevmode\nobreak\ t\geq 0,$ and the equation (17) can be transformed into [131]:

$${}\frac{\partial v^{s}}{\partial t}=\frac{\partial}{\partial x}\bigg{(}D(t-s)\frac{\partial v^{s}}{\partial x}+Bv^{s}\bigg{)}+d(t-s)v^{s},\leavevmode\nobreak\ \leavevmode\nobreak\ s\geq 0,\leavevmode\nobreak\ \leavevmode\nobreak\ s\leq t\leq s+\tau,$$

(18)

where $v^{s}(x,t)=u(x,t-s,t)$.

Nonlocal problems and conservation laws are particularly relevant in modelling biological transport and reaction-diffusion processes. They describe the movement of substances or entities within biological systems and the interactions between different components. Based on the mass conservation law, many mathematical models have been studied in cellular populations [132]. For example, if $u(\mathbf{x},t)$ and $J(\mathbf{x},t)$ are the density and flux of a population, then the conservation equation is given by

$${}u_{t}(\mathbf{x},t)=-\nabla\cdot J(\mathbf{x},t).$$

(19)

An appropriate choice of flux helps in capturing the different biological phenomena. As an example, the flux $J$ can be divided into the random motion ($J_{d}$) and the adhesion ($J_{a}$), i.e.,

$$J(\mathbf{x},t)=J_{d}(\mathbf{x},t)+J_{a}(\mathbf{x},t).$$

We assume Fick’s law for the diffusive flux, i.e., $J_{d}=-d_{1}\nabla u$. On the other hand, following Armstrong et al. [133], the adhesion flux is proportional to the population density and adhesion force. It is inversely proportional to the cell size, so

$$J_{a}=\frac{c}{L}u(\mathbf{x},t)F(\mathbf{x},t),$$

where $c$ is the proportionality constant, $F$ is the adhesion force, and $L$ is the cell size. Using Stoke’s law, the total adhesion force is given by

$$F(\mathbf{x},t)=\int_{\mathbb{B}^{n}(R)}g(u(\mathbf{x}+r,t))P(r)dr,$$

where $g(\cdot)$ and $P(\cdot)$ represent the nature of the force and its direction, respectively. Taking all these factors in the conservation equation (19), we obtain

$${}u_{t}(\mathbf{x},t)=d_{1}\Delta u(\mathbf{x},t)-\nabla\cdot(u(\mathbf{x},t)(\Phi\ast u)(\mathbf{x},t)),$$

(20)

where

$$(\phi\ast u)(\mathbf{x},t)=\frac{c}{L}\int_{\mathbb{B}^{n}(R)}g(u(\mathbf{x}+r,t))P(r)dr.$$

In [134], authors have considered the nonlocal model (20) with an additional term (logistic growth) in the reaction kinetics. They derived sufficient conditions for the boundedness of the solutions, but these conditions do not satisfy numerically for some parameter values, and hence, there are still some mathematical challenges.

In traditional variational problems, the goal is to find an extremum of a function that depends on local derivatives of a function. Nonlocal variational problems are mathematical optimization problems that involve nonlocal operators in their formulations. These formulations have been used to derive and study nonlocal PDEs in diverse physical and biological systems. The general form of a nonlocal variational problem can be expressed as follows:

$${}\mbox{Min or Max}\leavevmode\nobreak\ J(u)=\int_{\Omega}L(\mathbf{x},\mathbf{y},u(\mathbf{x}),u(\mathbf{y}))d\mathbf{x}d\mathbf{y}$$

(21)

subject to boundary conditions and constraints on the function $u(\mathbf{x})$ over the domain $\Omega$. Here, the function $L(\mathbf{x},\mathbf{y},u(\mathbf{x}),u(\mathbf{y}))$ is the nonlocal kernel that characterizes the interaction between the points $\mathbf{x}$ and $\mathbf{y}$ in the domain $\Omega$, and $u(\mathbf{x})$ represents the unknown function being optimized.

Nonlocal mathematical model also arises from crop-raiding of large-bodied mammals living in the biodiversity-rich tropics [135]. Let $\Omega$ be a spatial domain with $\overline{\Omega}_{0}\subset\Omega$ and $u(\mathbf{x},t)$ be the population density of the mammal at position $\mathbf{x}$ and time $t$. Here, we assume that the forest is safe for the mammal species but is of poor resources, while the farm has rich resources but is dangerous. Furthermore, mammals cannot survive if they only stay in the forest and do not attempt to go out to search for food; they cannot produce offspring on the farm because it is not a safe place for reproduction. In this case, the model takes the form [135]:

$${}u_{t}=d\Delta u+\gamma v(\mathbf{x})u\int_{\Omega\smallsetminus\Omega_{0}}u(\mathbf{y},t)d\mathbf{y}-u^{p},\leavevmode\nobreak\ \leavevmode\nobreak\ \mbox{in}\leavevmode\nobreak\ \Omega\times(0,T),$$

(22)

with Dirichlet boundary condition on $\partial\Omega$ and positive initial condition. Here, $\gamma>0$ and $v(\mathbf{x})$ is the characteristic function of $\Omega_{0}$, which is defined by

$$v(\mathbf{x})=\begin{cases}1\leavevmode\nobreak\ \leavevmode\nobreak\ \mbox{if}\leavevmode\nobreak\ \mathbf{x}\in\Omega_{0},\\ 0\leavevmode\nobreak\ \leavevmode\nobreak\ \mbox{if}\leavevmode\nobreak\ \mathbf{x}\in\Omega\smallsetminus\Omega_{0}.\end{cases}$$

The parameter $d$ is the diffusivity, and $p\geq 2$ takes into account the over-crowding effects of the population in $\Omega$. The non-negative steady-states for the model (22) are the nonnegative solutions of the semilinear nonlocal elliptic problem

$$-d\Delta u=\gamma v(\mathbf{x})u\int_{\Omega\smallsetminus\Omega_{0}}u(\mathbf{y},t)d\mathbf{y}-u^{p},\leavevmode\nobreak\ \leavevmode\nobreak\ \mbox{in}\leavevmode\nobreak\ \Omega\times(0,T).$$

A similar type of Kirchhoff-Carrier type equation boundary value problems studied in the literature [136]:

$$-(a+b\alpha(||u||,|u|_{\gamma}))\Delta u=\lambda\bigg{(}a_{1}+b_{1}\bigg{(}\int_{\Omega_{0}}c(x)|u|^{q_{1}}\bigg{)}^{\gamma_{1}}\bigg{)}u^{q_{2}}+f(u),$$

where $|u|_{\gamma}=(\int_{\Omega}|u|^{\gamma})^{1/\gamma}$, $||u||=(\int_{\Omega}|\nabla u|^{2})^{1/2}$, $\alpha$ is a non-negative function, and $c$ is a non-negative continuous function in $\Omega$. The conditions on $f$ and the other parameters are given in [136]. A weak solution to this problem is a function $u$ in a Hilbert space $H_{0}^{1}(\Omega)$ such that

$$(a+b\alpha(||u||,|u|_{\gamma}))\int_{\Omega}\nabla u\cdot\nabla h=\lambda\bigg{(}a_{1}+b_{1}\bigg{(}\int_{\Omega_{0}}c(x)|u|^{q_{1}}\bigg{)}^{\gamma_{1}}\bigg{)}\int_{\Omega}u^{q_{2}}h+\int_{\Omega}f(u)h,\leavevmode\nobreak\ \leavevmode\nobreak\ \forall h\in H_{0}^{1}(\Omega).$$

Furthermore, the variational structure $J_{\lambda}(u)$ has derivative operator given by

$$J_{\lambda}^{\prime}(u)h=\bigg{(}\int_{\Omega_{0}}c(x)|u|^{q_{1}}\bigg{)}^{\gamma_{1}}\int_{\Omega}u^{q_{2}}h,\leavevmode\nobreak\ \leavevmode\nobreak\ \forall u,h\in H_{0}^{1}(\Omega).$$

Structure-preserving characteristics and variational integrators are fundamental notions in numerical analysis, particularly in the simulation of dynamical systems [137]. These approaches are intended to retain fundamental aspects of the underlying physical system, such as energy conservation, ensuring that the numerical solution closely resembles the genuine system’s behaviour over long time periods. These methodologies are used in a variety of scientific areas, ranging from celestial mechanics to molecular dynamics simulations.

Sometimes, the complete information of a system’s states may be inaccessible, but only partial information about the states is observable due to limitations in experimental techniques. In this case, observability is a critical aspect of partially observed dynamics, highlighting which parts of a system are accessible to measurement [138]. The reduction techniques focus on accounting for hidden variables and incorporating partial observations into the system’s modelling. Stochastic models and Bayesian inference are employed to handle uncertainties associated with partially observed dynamics, and these frameworks provide probabilistic representations of hidden states and update beliefs based on observed data.

Pattern analysis was essential to human survival and advancement in the past. In order to make sense of the world, ancient civilizations studied patterns in astronomy, agriculture, social behaviour, and nature. Seasonal trends and natural cycles were intimately related to agricultural operations. To keep the soil fertile and minimize depletion, farmers would rotate the crops they grew in a given field. Numerous ancient societies, like the Greeks, Egyptians, and Mesopotamians, depended on the yearly inundation of rivers like the Nile and Tigris-Euphrates to replenish the soil with nutrients and used astronomical patterns as a guide. Although the first agricultural communities probably lacked the sophisticated maps we have today, they did manage and plan their agricultural activities using simple spatial organizing techniques. For instance, the historic Indus Valley civilization metropolis of Mohenjo-Daro which dates to approximately 2600 BCE. Cuneiform tablets uncovered in ancient Mesopotamia, one of the first centres of civilization, date to approximately 2300 BCE and show primitive maps of land property.

History demonstrates that the time of agricultural activities was greatly influenced by the locations of the sun, moon, and stars. A Greek mathematician and astronomer, Hipparchus (190–120 BCE), is frequently acknowledged for having made some of the first attempts at star mapping. By mapping the stars using a coordinate system, he popularized the ideas of celestial longitude and latitude. Hipparchus’s suggestion of a method to measure the apparent brightness of stars is one of his most important contributions. Babylonian mathematicians also investigated natural numerical patterns. Pythagorean triples are listed on the Plimpton 322 clay tablet from ancient Babylonia, which dates to approximately 1800 BCE and illustrates mathematical patterns in ancient times. In addition, traditional medical professionals and healers examined trends in disease and symptomology. They helped to shape the early forms of herbal medicine and other medical practices by observing the impact of specific plants or ceremonies on health. Ancient civilizations didn’t have the advanced technology that we use today, but they did develop their own methodologies to understand the world. Their models may not have been explicitly nonlocal in the way we think about it today, but they certainly reflected an understanding of interconnectedness and patterns in the world.

Nonlocal interactions in ecological systems can have a significant impact on spatiotemporal pattern formation, e.g., the development of spatial patterns in plants. Nonlocal interactions between plants, such as resource competition or facilitation, can result in the creation of regular patterns such as stripes, spots, or waves. These patterns frequently emerge as a result of self-organization in which interactions between individuals on a broader scale give rise to organized patterns on a landscape level. Nonlocal interactions in predator-prey relationships can result in the development of predator fronts, in which predators pursue prey concentrations over longer distances, altering the spatial dynamics of both populations. The nonlocal predator-prey interactions are generally two types: intra- and inter-specific. First, we consider a general nonlocal prey-predator model with nonlocal interaction in the intra-specific prey competition in a two-species population as

	$\displaystyle\frac{\partial u}{\partial t}$	$\displaystyle=d_{1}\frac{\partial^{2}u}{\partial x^{2}}+\alpha u\bigg{(}1-\frac{\phi\ast u}{\kappa}\bigg{)}-f(u)v,$		(23)
	$\displaystyle\frac{\partial v}{\partial t}$	$\displaystyle=d_{2}\frac{\partial^{2}v}{\partial x^{2}}+(ef(u)-g(v))v,$

subjected to non-negative initial conditions and appropriate boundary conditions of prey ($u$) and predator ($v$) populations. The parameters $\alpha$ and $\kappa$ are intrinsic growth rate and environmental carrying capacity, respectively, and $e$ $(0<e\leq 1)$ is the conversion efficiency from prey biomass into predator biomass. If $e=0$, then there is no interaction between the species. For simplicity, we have considered a one-dimensional spatial domain $\Omega\subset\mathbb{R}$, and $d_{1}$ and $d_{2}$ are the self-diffusion coefficients. As defined earlier, $\phi\ast u$ is the convolution term. The kernel function $\phi$ satisfies all the assumptions mentioned earlier, i.e., normalized, and has compact support in $\mathbb{R}$. For an unbiased movement of the species in both spatial directions, we assume $\phi$ as an even function. Still, sometimes, species do not have the same area for nonlocal consumption of resources in both spatial directions. This generally happens for a bounded spatial domain, e.g., studying a nonlocal interaction for the populations living in a river [46, 139, 140]. The function $f(u)$ is called a functional response, and it satisfies three conditions [141, 142]: (i) $f(0)=0$, (ii) $f$ is an increasing function and (iii) there exists a finite $M>0$ such that $\lim_{u\rightarrow\infty}f(u)=M$. Many biologically relevant functional responses are available in the literature, e.g., Holling type II and III, Ivlev function, etc. [141, 142, 143, 144]. The function $g(v)$ represents the per capita death rate of the predator population [145]. In [76], authors have studied a particular form of the nonlocal model (23). Nevertheless, the local model corresponding to that nonlocal model does not produce any nonhomogeneous stationary pattern, but the nonlocal model does.

Researchers have also studied the pattern formation for nonlocal models with more than two species populations. Generally, the three components of the local or nonlocal model in ecology add another feature to the two-species model. It allows us to study the dynamics of the prey-predator model in the presence of an additional species. These three components of the model may be (a) food chain; (b) two predators–one prey; (c) one predator–two prey; (d) food chain with omnivory; and (e) food chain with cycle [146]. Autry et al. [147] have considered a three-species nonlocal RD model in a food chain system with ratio-dependent functional responses as:

	$\displaystyle\frac{\partial u}{\partial t}$	$\displaystyle=d_{1}\frac{\partial^{2}u}{\partial x^{2}}+u\bigg{(}1-\phi\ast u-\frac{a_{1}v}{u+v}\bigg{)},$		(24)
	$\displaystyle\frac{\partial v}{\partial t}$	$\displaystyle=d_{2}\frac{\partial^{2}v}{\partial x^{2}}+v\bigg{(}-\mu_{1}+\frac{m_{1}u}{u+v}-\frac{a_{2}w}{v+\psi\ast w}\bigg{)},$
	$\displaystyle\frac{\partial w}{\partial t}$	$\displaystyle=d_{3}\frac{\partial^{2}w}{\partial x^{2}}+w\bigg{(}-\mu_{2}+\frac{m_{2}v}{v+\psi\ast w}\bigg{)},$

where $u$, $v$, and $w$ represent the concentration of prey, predator, and super-predator, respectively. All the convolution terms are defined as before, and $\phi$ and $\psi$ are the kernel functions. Here, the prey and super-predator interact directly with the intermediate predator and only interact indirectly with the intermediate predator. Authors in [147] have studied two aspects of the nonlocality: C-type and P-type. The $C$-type nonlocality occurs when the prey species competes with itself. This type of nonlocal interaction usually occurs due to the overlapping ranges of the species, e.g., extensive root systems or a shared water table in plant species. On the other hand, P-type stands for pest type, which arises when one species nonlocally interacts with the other species, and this nonlocal interaction only affects the encounter rate, which appears in the functional response term. In addition, one can consider the same type of nonlocal term in the functional response for the first and second equations to study the nonlocal effect of high-mobile intermediate predator species $v$ on the crops or prey $u$.

In population biology, age structure models and pattern formation provide insights into the dynamics of populations and their structural changes over time. Age structure refers to the distribution of individuals in a population across different age groups, and understanding its patterns is crucial for demographic analysis and prediction. McKendrick and Von Foerster developed this concept [148, 149]. Suppose $n(a,t)$ is the population at age $a$ and time $t$. Following McKendrick - Von Foerster, the population $n(a,t)$ satisfies the following equation [150]:

$$\frac{\partial n}{\partial t}+\frac{\partial n}{\partial a}+\mu(a)n(a,t)=0,$$

with the initial condition $n(a,0)=f(a)$ and the birth boundary condition:

$$n(0,t)=\int_{0}^{\infty}b(a)n(a,t)da,$$

where $b(a)$ is the birth rate of the individuals which depends on the age of the adult population. Here, $\mu(a)$ is the death rate of individuals. These types of age structure models often use population pyramids to represent the distribution of age groups in a population visually.

In pattern formation, morphogenesis is a process by which spatial patterns and structures arise during the development of organisms [151]. The creation of the various and complex forms found in organisms depends on this complex process. It involves a cascade of molecular and cellular processes that coordinate the orientation, polarity, and assembly of cells to generate distinct morphologies and configurations. Examples include (i) the development of an embryo into a fully formed organism, (ii) the differentiation process by which cells in developing organisms take on certain destinies and functions, and (iii) precise control of gene expression. The multidisciplinary study of morphogenesis aims to comprehend the fundamental processes that give rise to the astounding diversity of biological forms. It draws on the disciplines of genetics, developmental biology, and biophysics, among others.

Spatiotemporal patterns in epidemiology reveal the dynamic interplay of place, time, and disease distribution within populations. Environmental, social, and biological factors all have an impact on these patterns [152, 153]. Epidemics can follow seasonal trends, exhibit cyclic behaviour, or demonstrate irregular patterns driven by a variety of factors, such as host susceptibility and environmental conditions, which are referred to as temporal patterns. However, spatio-temporal models deal with many spatial heterogeneity factors, such as higher or lower incidence rates in certain areas [154].

Nonlocal models are essential for describing developing patterns in a variety of biological systems, especially when it comes to cell communication, cell adhesion, and bacterial aggregation. Cell adhesion models take into account distant adhesive forces, enabling cells to recognize and respond to the spatial distribution of adhesion molecules [133]. It can capture the collective behaviour of bacteria by addressing long-range interactions, such as quorum sensing, where bacterial cells communicate and coordinate their activities across longer distances [155].

Nonlocal models derived from quantum confinement are intriguing in characterizing pattern generation at the nanoscale and have several applications in materials research, quantum mechanics, and device manufacturing. Quantum confinement is the phenomenon in which the electronic and optical characteristics of materials are significantly affected by their size, particularly when limited to dimensions equivalent to electrons’ characteristic length scale [156]. The nonlocal models derived from quantum mechanics are frequently used to characterize the electrical and optical behaviour of quantum dots, nanowires, and other nanostructures [157, 158]. The discrete energy levels and wave functions of confined electrons produce distinct patterns in their characteristics. Quantum dot solar cells, light-emitting diodes, and quantum dot transistors are some of the uses [159, 160]. Nonlocal models can assist in explaining the development of energy bands and electronic patterns in quantum dot arrays, which influence their optical and electrical characteristics [161]. Furthermore, quantum confinement can cause nonlocal electron interactions and transport behaviours, affecting the operation of nanoelectronic devices such as sensors, catalysis, and quantum computing.

An interesting perspective on the emergence of complex patterns in many physical and biological systems may be gained from nonlocal models that describe pattern development resulting from symmetry-breaking processes [162, 163]. Symmetry-breaking is a basic idea in which complex patterns are formed when a system changes from a symmetrical to an asymmetrical state. These occurrences, prompted by disturbances or outside factors, start to establish patterns. Fluid dynamics, chemical reactions, and biological processes are examples of systems where this may happen [164, 165, 166, 167, 168].