A Physical Phenomenon for the Fractional Nonlinear Mixed Integro-Differential Equation with Local and Nonlocal Conditions Using a Toeplitz Matrix Technique with a Genetic Application

1. Introduction

Fractional differential equations have garnered considerable interest in recent years owing to their capacity to simulate intricate systems characterized by memory and hereditary attributes. Unlike classical differential equations, fractional equations involve derivatives and integrals of non-integer order, which allow for more accurate descriptions of various physical and engineering phenomena [1]. In particular, fractional integral equations with nonlocal boundary conditions have emerged as an essential area of study, with applications in viscoelasticity, anomalous diffusion, and biological systems.

Existence and uniqueness of solutions of fractional integral and differential equations with nonlocal constraints have been studied in several works. Ahmad and Luca [2] focused on a system of fractional differential equations with coupled nonlocal boundary conditions and derived existence results using fixed-point theorems. Similarly, Gou and Li [3] examined impulsive Hilfer fractional evolution equations and established theoretical results for nonlocal problems.

From a numerical perspective, various methods have been developed to approximate solutions to fractional problems. Jiang et al. [4] proposed a modified minimum residual method for nonlinear fractional equations with nonlocal boundary conditions, demonstrating its effectiveness in solving complex systems. Sultanov et al. [5] introduced a numerical method for subdiffusion differential equations with nonlocal boundary conditions, further contributing to computational approaches in this field.

There has also been some interest in stability analysis of fractional systems with nonlocal criteria. A nonlinear Hadamard fractional coupling Laplacian system with symmetric periodic boundary conditions was studied in stability and numerical simulation by Lv et al. [6]. In the meanwhile, Ye and Qu [7] examined non-local solutions to fractional non-instantaneous impulse integro-differential equations that involve mild solutions.

Additionally, functional fractional equations with infinite delay and nonlocal integral boundary conditions have been studied by Alghanmi and Alqurayqiri [8], extending the applicability of fractional models to dynamic systems with memory effects. Further, Mosal et al. [9] discussed the solution of mixed IE with nonlocal conditions. Raad [10] discussed the solution of local and non-local phase- lag integro partial differential solution.

From an alternative perspective, evaluate and contrast a variety of numerical solution and research strategies, with the fundamental concept of the equation serving as the determining factor in the selection of the most effective and resilient approach. The current investigation will evaluate the research that has implemented these methodologies. The quadrature Nyström technique for fractional nonlinear mixed integro-differential equations was devised by Jan et al. [11] and its numerical efficacy was proved. The quadratic phase-lag integral problems were theoretically and computationally examined by Abdou et al. [12]. Rostami and Maleknejad [13] employed hybrid numerical methods to resolve two-dimensional nonlinear mixed partial integro-differential equations. They employed modified block-pulse functions and Bernoulli polynomials. The Galerkin method was implemented by Sayed et al. [14] to identify solutions to hyperbolic telegraph equations in one and two dimensions. The Legendre spectral collocation method was investigated by Ezz-Eldien and Alalyani [15] for the two-dimensional Volterra-Fredholm integral problem. The numerical solutions of the nonlinear ordinary and fractional Newell-Whitehead-Segel equations were investigated by Yassin et al. [16]. Abd-Elhameed et al. [17] have employed the seventh kind of Chebyshev polynomials approach to resolve fractional delay differential equations.

In this paper, we aim to contribute to the theoretical and numerical analysis of FrNMIo-DE with nonlocal conditions. Local conditions are initial conditions in the absence of time, while in nonlocal conditions, the time condition is taken into account before the initial condition. The concept of non-local condition in mathematical biology is the effect of genetic state on the bodies being studied, and at the same time, the time delay is considered, meaning that the non-local condition is observed for its effects in genetic engineering, whether in humans or agricultural materials, by taking into account the previous time before its emergence.

We explore existence and uniqueness results, provide stability analysis, and propose numerical algorithms for solving such equations. Our approach builds upon previous works while introducing novel techniques to enhance accuracy and computational efficiency. In the second section, the fractional-time differential-integral equation was reviewed, taking into account nonlocal conditions. The resulting integral equation was then calculated using the fundamental rules of fractional integrals. Finally, the conditions necessary for obtaining a unique solution were established. In the third section, the existence of a unique solution was investigated using the fixed-point theorem. In the fourth section, the authors addressed the convergence of the solution and proved it to be uniform. The fifth section employed the position-time separation method. The sixth section explored the compact matrix method for solving the problem numerically, outlining the fundamental reasons for its preference over other methods. The seventh section examined the unity of the nonlinear algebraic system obtained using the compact matrix method. The eighth section investigated the convergence of this nonlinear algebraic system. Finally, the ninth section analyzed the resulting error and its convergence rate. In section ten, some applications were solved numerically when the nucleus had certain special shapes, such as logarithmic nuclei and Carlman nuclei, while studying the effect of local and non-local conditions between equations. Section 11 describes a physical Justification: Local and Nonlocal Conditions in Genetic Gene Regulation. Section 12. A general conclusion of the numerical results presented in this paper provide a complete and systematic proof of the accuracy, stability and convergence of the proposed separation method coupled with the Toplitz matrix method for the numerical solution of fractional integral-differential equations with both local and non-local initial conditions.

2. Time Fractional and Nonlinear Integro-Differential Equation

Consider, in the time fractional science, the following integro-differential equation

$$\left({\mu }_0 {\displaystyle \frac{{\partial }^{\alpha +1}}{\partial t^{\alpha +1}}}+{\mu }_1{\displaystyle \frac{{\partial }^{\alpha }}{\partial t^{\alpha }}}+{\mu }_2\right)\phi \left(x,t\right)=g\left(x,t\right)+\mathrm{br}-\mathrm{to}-\mathrm{break} \lambda F\left(t\right){\int }_{-1}^{1}k\left(\left|x-y\right|\right)\gamma \left(y,\tau ,\phi \left(y,\tau \right)\right)dy,\left(0<\mathsf{\alpha }<1\right).$$

(1)

Under the nonlocal conditions

$$\phi \left(x,0\right)=v_1\left(x\right)+w\left(x,t\right),{\displaystyle \frac{\partial \phi \left(x,0\right)}{\partial t}}=v_2\left(x\right)+w_{\mathrm{t}}\left(x,t\right).$$

(2)

In this context, g(x,t) is a defined function within the space ${\mathrm{L}}_2[-\mathrm{1,1}]\times \mathrm{C}[0,\mathrm{T}]$, where $\mathrm{T}<1$, serving as the free term of the issue. The constant $\mathsf{\lambda }$ is contingent upon the material type and holds several physical significances in applied mathematics. The function $\mathrm{F}\left(\mathrm{t}\right)$serves as a smooth kernel in time, whereas $k\left(\right|x-y\left|\right)$is a single kernel in position, represented in logarithmic form and as a Carleman function. $\mathsf{\gamma }(\mathrm{x},\mathrm{t},\mathsf{\phi }(\mathrm{x},\mathrm{t}\left)\right)$is a specified nonlinear function of the unknown function $\phi (x,t). v_1\left(x\right)$ and $v_2\left(x\right)$ are two specified beginning position functions. $\left\{w\right(x,t), w_t(x,t\left)\right\}$ are temporal and spatial historical functions. By integrating, equation (1) becomes

$${\mathsf{\mu }}_0{\displaystyle \frac{{\partial }^{\alpha }}{\partial t^{\alpha }}}\left(\phi \left(x,t\right)-{\displaystyle \frac{\partial }{\partial t}}\phi \left(x,0\right)\right)+{\mathsf{\mu }}_1{I}^{1-\alpha }\phi \left(x,t\right)+{\mathsf{\mu }}_2{\int }_0^t\phi \left(x,\tau \right)d\tau$$

$$={\int }_0^{t}g\left(x,\tau \right)d\tau +\lambda {\int }_0^t{\int }_{-1}^{1}F\left(\tau \right)k\left(\left|x-y\right|\right)\gamma \left(y,\tau ,\phi \left(y,\tau \right)\right)dyd\tau .$$

(3)

Using the basic formula of Caputo-fractional integral

$$I_a^{\alpha }f\left(x\right)={\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\int }_a^x{\left(x-t\right)}^{\alpha -1}f\left(t\right)dt,$$

(4)

$${\int }_0^t{\int }_0^{{\tau }_{n-1}}\dots {\int }_0^{{\tau }_2}{\int }_0^{{\tau }_1}f\left(\tau \right)d\tau d{\tau }_1\dots d{\tau }_{n-2}d{\tau }_{n-1}={\displaystyle \frac{1}{\mathsf{\Gamma }\left(n\right)}}{\int }_0^t{\left(t-\tau \right)}^{n-1}f\left(\tau \right)d\tau ,$$

(5)

equation (3) tends to

$${\mu }_0\phi \left(x,t\right)+{\int }_0^t\left({\mu }_1+{\displaystyle \frac{{\mu }_2}{\mathsf{\Gamma }\left(\alpha \right)}}{\left(t-\tau \right)}^{\alpha }\right)\phi \left(x,\tau \right)d\tau =G\left(x,t\right)+{\displaystyle \frac{\mathsf{\lambda }}{\mathsf{\Gamma }\left(\alpha \right)}}{\int }_0^t{\int }_{-1}^1{\left(t-\tau \right)}^{\alpha }F\left(\tau \right)k\left(\left|x-y\right|\right)\gamma \left(y,\tau ,\phi \left(y,\tau \right)\right)dyd\tau ,$$

$$G\left(x,t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\int }_0^t{\left(t-\tau \right)}^{\alpha }g\left(x,\tau \right)d\tau -{\displaystyle \frac{{\mu }_0{\mathrm{t}}^{\alpha }}{\alpha \Gamma \left(\alpha \right)}}\left({\mathrm{v}}_1\left(\mathrm{x}\right)+{\mathrm{v}}_2\left(\mathrm{x}\right)\right)+{\displaystyle \frac{{\mu }_0}{\Gamma \left(\alpha \right)}}{\int }_0^t{\left(t-\tau \right)}^{\alpha -1}\left(\mathrm{w}\left(\mathrm{x}, \tau \right)+{\mathrm{w}}_t\left(\mathrm{x}, \tau \right)\right)d\tau .$$

(6)

The following conditions are assumed in order to ensure the existence of a unique solution to the problem of Equation (1) or its equivalent in Equation (6):

i) The norm of the unknown function$\phi (x,t)$and its derivatives is defined as, and it is located in the space$L_2 [-\mathrm{1,1}]\times C[0,T].$

$$\Vert \phi \left(x,t\right)\Vert =\underset{0\le t\le T}{\max }{\int }_0^t{\left\{{\displaystyle \underset{-1}{\overset{1}{\int }}}{\phi }^2\left(x,\tau \right)dx\right\}}^{1/2}d\tau .$$

ii) For the constant$A$, the known function $\gamma \left(x,t,\phi \right)$ satisfies the following conditions:

ii-a)$\Vert \gamma \left(x,t,\phi \right)\Vert \le A\Vert \phi \left(x,t\right)\Vert .$ii-b) $\Vert \gamma \left(x,t,{\phi }_1\right)-\gamma \left(x,t,{\phi }_2\right)\Vert \le A\Vert {\phi }_1\left(x,t\right)-{\phi }_2\left(x,t\right)\Vert .$iii) The given function $g(x,t)$ satisfies

$$\Vert g\left(x,t\right)\Vert =\underset{0\le t\le T}{\max }{\int }_0^t{\left\{{\displaystyle \underset{-1}{\overset{1 }{\int }}}{g}^2\left(x,\tau \right)dx\right\}}^{1/2}d\tau \le B,B\mathrm{i}\mathrm{s}\mathrm{a}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t}.$$

iv) For the constants $P_i, i=\left\{\mathrm{1,2}\right\}$, the two functions $v_i\left(x\right)$ fulfil the following conditions:

iv-a) $\left|v_1\left(x\right)\right|\le P_1$ and $|v_2 (x\left)\right|\le P_2.$For the constants $Q_i , i=\left\{1, 2\right\},$ $w\left(x,t\right)$ and $w_t\left(x,t\right)$ satisfy

iv-b) $\Vert w\left(x,t\right)\Vert =\underset{\begin{array}{c}0\le t\le T,\\ \end{array}}{\mathrm{m}\mathrm{a}\mathrm{x}} {\int }_0^t{\int }_{-1}^1{{\{w\left(x,\tau \right)}^{2}dx\}}^{{\displaystyle \frac{1}{2}}}d\tau \le Q_1$,

iv-c) $\Vert w_{\mathrm{t}}\left(x,t\right)\Vert =\underset{\begin{array}{c}0\le t\le T,\\ -1\le x\le 1\end{array}}{\mathrm{m}\mathrm{a}\mathrm{x}} \left|w\left(x,t\right)-w\left(x,0\right)\right|\le Q_2$v) The position kernel in the space $L_2[-\mathrm{1,1}]$ satisfies

$$\Vert k\left(\left|x-y\right|\right)\Vert ={\left\{{\int }_{-1}^1{\int }_{-1}^1{k}^2\left(\left|x-y\right|\right)dxdy\right\}}^{1/2}=C.$$

vi) The continuous function $F\left(t\right)$ in time satisfies

$$\underset{0\le t\le T}{\mathrm{m}\mathrm{a}\mathrm{x}}{\int }_0^{t}F\left(\tau \right)d\tau =E$$

3. Existence and Uniqueness

Equation (6) can be expressed in the following integral operator form to demonstrate its existence and uniqueness.

$$\chi \phi \left(x,t\right)={\chi }_1\phi \left(x,t\right)-{\chi }_2\phi \left(x,t\right)+{\displaystyle \frac{1}{{\mu }_0\Gamma \left(\alpha \right)}}{\int }_0^t{\left(t-\tau \right)}^{\alpha }g\left(x,\tau \right)d\tau -{\displaystyle \frac{{\mu }_0{\mathrm{t}}^{\alpha }}{\Gamma \left(\alpha +1\right)}}\left({\mathrm{v}}_1\left(\mathrm{x}\right)+{\mathrm{v}}_2\left(\mathrm{x}\right)\right)+{\displaystyle \frac{{\mu }_0}{\Gamma \left(\alpha \right)}}{\int }_0^t{\left(t-\tau \right)}^{\alpha -1}\left(\mathrm{w}\left(\mathrm{x}, \tau \right)+\left(\mathrm{x}, \tau \right)\right)d\tau ,$$

$${\chi }_1\phi \left(x,t\right)={\displaystyle \frac{\mathsf{\lambda }}{{\mu }_0\Gamma \left(\alpha \right)}}{\int }_0^t{\int }_{-1}^1{\left(t-\tau \right)}^{\alpha }F\left(\tau \right)k\left(\left|x-y\right|\right)\gamma \left(y,\tau ,\phi \left(y,\tau \right)\right)dyd\tau ,{\chi }_2\phi \left(x,t\right)={\displaystyle \frac{1}{{\mu }_0}}{\int }_0^t\left({\mu }_1+{\displaystyle \frac{{\mu }_2}{\mathsf{\Gamma }\left(\alpha \right)}}{\left(t-\tau \right)}^{\alpha }\right)\phi \left(x,\tau \right)d\tau .$$

(7)

Theorem 1. 

The solution of Equation (7) exists and unique under the condition

$$\left|\lambda \right|T^{\alpha +1}ECA+\left|{\mu }_1\right|\mathrm{T}\mathsf{\Gamma }\left(\alpha +1\right)+\left|{\mu }_2\right|T^{\alpha +1}<\left|{\mu }_0\right|\mathsf{\Gamma }\left(\alpha +1\right).$$

(8)

The following lemmas must be proved to satisfy the above theorem.

The operator χ transforms the space${\mathrm{L}}_2[-\mathrm{1,1}]\times \mathrm{C}[0,\mathrm{T}]$

onto itself based on criteria (i)-(vi).

Proof: From Equation (7), we get

$$\Vert \chi \phi \left(x,t\right)\Vert \le \Vert {\chi }_1\phi \left(x,t\right)\Vert +\Vert {\chi }_2\phi \left(x,t\right)\Vert +{\displaystyle \frac{1}{\left|{\mu }_0\right|\Gamma \left(\alpha \right)}}\left|{\int }_0^t{\left(t-\tau \right)}^{\alpha }g\left(x,\tau \right)d\tau \right|+{\displaystyle \frac{t^{\alpha }}{\Gamma \left(\alpha +1\right)}}\left(\left|v_1\left(x\right)\right|+\left|v_2\left(x\right)\right|\right)+ {\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}\left|{\int }_0^t{\left(t-\tau \right)}^{\alpha -1}\left(\mathrm{w}\left(\mathrm{x}, \tau \right)+{\mathrm{w}}_t\left(\mathrm{x}, \tau \right)\right)d\tau \right|.$$

(9)

Using conditions (i)-(v) and Cauchy-Schwartz inequality, we have

$$\Vert \chi \phi \left(x,t\right)\Vert \le \delta \Vert \phi \left(x,t\right)\Vert +{\displaystyle \frac{\mathrm{B}{T}^{\alpha +1}}{\left|{\mu }_0\right|\mathsf{\Gamma }\left(\alpha +1\right)}}+{\displaystyle \frac{{\mathrm{T}}^{\alpha }}{\mathsf{\Gamma }\left(\alpha +1\right)}}\left(P_1+{P_2+Q}_1+Q_2\right),$$

(10)

Where

$$\delta ={\displaystyle \frac{\left|\lambda \right|T^{\alpha +1}ECA+\left|{\mu }_1\right|\mathrm{T}\mathsf{\Gamma }\left(\alpha +1\right)+\left|{\mu }_2\right|T^{\alpha +1}}{\left|{\mu }_0\right|\mathsf{\Gamma }\left(\alpha +1\right)}}.$$

(11)

It is obvious that the operator $\chi$ maps the ball $B_r \in L_2\left[-1, 1\right]\times C\left[0, T\right]$Š into itself where

$r={\displaystyle \frac{\sigma }{1-\delta }}$ , σ =${\displaystyle \frac{\mathrm{B}{T}^{\alpha +1}}{\left|{\mu }_0\right|\mathsf{\Gamma }\left(\alpha +1\right)}}+{\displaystyle \frac{{\mathrm{T}}^{\alpha }}{\mathsf{\Gamma }\left(\alpha +1\right)}}\left(P_1+{P_2+Q}_1+Q_2\right)$.

The inequality (10) involves the boundedness of the operator $\chi$ under the condition δ < 1.

Lemma 2 

. If requirements (i)–(vi) are fulfilled, then χ is a contraction operator in Banach space $L_2\left[-1, 1\right]\times C\left[0, T\right].$

Proof. Let two functions ${\phi }_1\left(x,t\right)$ and ${\phi }_2\left(x,t\right)$ be two solutions of (6), and then, the formula (7) leads to$the$ formula

$$\Vert \chi {\phi }_1\left(x,t\right)-\chi {\phi }_2\left(x,t\right)\Vert \le \Vert {\chi }_1\left({\phi }_1\left(x,t\right)-{\phi }_2\left(x,t\right)\right)\Vert +\Vert {\chi }_2\left({\phi }_1\left(x,t\right)-{\phi }_2\left(x,t\right)\right)\Vert .$$

Using conditions (i)-(iv) and Cauchy Schwarz inequality, we deduce that

$$\Vert \chi {\phi }_1\left(x,t\right)-\chi {\phi }_2\left(x,t\right)\Vert \le \delta \Vert {\phi }_1\left(x,t\right)-{\phi }_2\left(x,t\right)\Vert$$

(12)

Therefore, for δ < 1, χ is a contraction operator in Equation (7). As a result, a unique solution in L_2 [[-1,1] exists for every $\mathrm{t}\in \mathrm{C}\left[0,\mathrm{T}\right],\mathrm{T}<1$, according to the Banach fixed point theorem.

4. Convergence of Solution

For this aim, take the straightforward iteration $\left\{{\phi }_1\left(x, t\right), {\phi }_2\left(x, t\right), ..., {\phi }_n\left(x, t\right), ... \right\}\subset \phi \left(x, t\right).$ Then, use Equation (7), to have

$${\mu }_0\left({\phi }_n\left(x,t\right)-{\phi }_{n-1}\left(x,t\right)\right)+{\int }_0^t\left({\mu }_1+{\displaystyle \frac{{\mu }_2}{\mathsf{\Gamma }\left(\alpha \right)}}{\left(t-\tau \right)}^{\alpha }\right)\left({\phi }_{n-1}\left(x,\tau \right)-{\phi }_{n-2}\left(x,\tau \right)\right)d\tau ={\displaystyle \frac{\lambda }{\Gamma \left(\alpha \right)}}{\int }_0^t{\int }_{-1}^1{\left(t-\tau \right)}^{\alpha }F\left(\tau \right)k\left(\left|x-y\right|\right)\left({\gamma }_n\left(y,\tau ,{\phi }_{n-1}\left(y,\tau \right)\right)-{\gamma }_n\left(y,\tau ,{\phi }_{n-2}\left(y,\tau \right)\right)\right)dyd\tau$$

(13)

Consider

$${\phi }_n\left(x,t\right)={\displaystyle \sum _{i=0}^n}{\xi }_i\left(x,t\right),$$

where

$${\xi }_n\left(x,t\right)={\phi }_n\left(x,t\right)-{\phi }_{n-1}\left(x,t\right), \left(n\ge 1\right).$$

and

$${\xi }_0\left(x,t\right)={\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_0^t{\left(t-\tau \right)}^{\alpha }g\left(x,\tau \right)d\tau +{\displaystyle \frac{{\mu }_0{t}^{\alpha }}{\alpha \Gamma \left(\alpha \right)}}\left(v_1\left(x\right)+v_2\left(x\right)\right) {\displaystyle \frac{{\mu }_0}{\Gamma \left(\alpha \right)}}{\int }_0^t{\left(t-\tau \right)}^{\alpha -1}\left(\mathrm{w}\left(\mathrm{x}, \tau \right)+{\mathrm{w}}_t\left(\mathrm{x}, \tau \right)\right)d\tau .\left(14\right)$$

(14)

Lemma 3. 

A sequence $\left\{{\mathsf{\phi }}_{\mathrm{n}}\left(\mathrm{x},\mathrm{t}\right)\right\}$

of Equation (14) is uniformly convergent under the condition $\delta <1.$

Proof: By applying Cauchy-Schwarz and using (14) in (13), we get

$$\left|{\mu }_0\right|\Vert {\xi }_n\left(x,t\right)\Vert \le \left|{\int }_0^t\left({\mu }_1+{\displaystyle \frac{{\mu }_2}{\mathsf{\Gamma }\left(\alpha \right)}}{\left(t-\tau \right)}^{\alpha }\right)d\tau \right|\Vert {\xi }_{n-1}\left(x,t\right)\Vert$$

(15)

$$+{\displaystyle \frac{\lambda A}{\Gamma \left(\alpha \right)}}\left|{\int }_0^t{\int }_{\Omega }{\left(t-\tau \right)}^{\alpha }F\left(\tau \right)k\left(\left|x-y\right|\right)dyd\tau \right|\Vert {\xi }_{n-1}\left(x,t\right)\Vert$$

Taking $n=1$, the above formula become

$$\Vert {\xi }_1\left(x,t\right)\Vert \le \delta \left[{\displaystyle \frac{\mathrm{B}{T}^{\alpha +1}}{\mathsf{\Gamma }\left(\alpha +1\right)}}+{\displaystyle \frac{\left|{\mu }_0\right|{\mathrm{T}}^{\alpha }}{\alpha \mathsf{\Gamma }\left(\alpha \right)}}\left(P_1+{P_2+Q}_1+Q_2\right)\right].$$

And, in general we have then

$$\Vert {\xi }_n\left(x,t\right)\Vert \le {\delta }^n\left[{\displaystyle \frac{\mathrm{B}{T}^{\alpha +1}}{\mathsf{\Gamma }\left(\alpha +1\right)}}+{\displaystyle \frac{\left|{\mu }_0\right|{\mathrm{T}}^{\alpha }}{\alpha \mathsf{\Gamma }\left(\alpha \right)}}\left(P_1+{P_2+Q}_1+Q_2\right)\right],\delta <1.$$

(16)

The sequence $\left\{{\xi }_n\left(x,t\right)\right\}$ is uniformly convergent by Equation (16). Additionally, it gives the sequence's ${\xi }_n\left(x,t\right)={\displaystyle \sum _{i=0}^n}{\xi }_i\left(x,t\right)$ convergent solution. As $n\to \infty , {\phi }_n\left(x,t\right)\to \phi \left(x,t\right)$, hence the solution $\phi \left(x,t\right)$ is uniformly convergent under the condition $\delta <1$. This demonstrates the lemma.

5. Separation of Variables Technique

In mathematical physics problems, the undiscovered possible function, related to space and time, is the object of researchers' attention. There are several ways to find the unknown function. One such method is time division, which reduces the mixed integral issue to an algebraic set of integral equations.

A sophisticated mathematical tool, the separation technique approach is used by certain researchers to translate FrNIE with time-dependent into a class of integral equations with time-only coefficients. By decoupling the integral equations from the time-dependent component, this method streamlines the problem and makes it easier to solve. It is possible to see the known and unknown functions separated in the form as

$$\phi \left(x,t\right)=N\left(t\right)\psi \left(x\right),f\left(x,t\right)=g\left(x\right)M\left(t\right),\gamma \left(x,t,\phi \left(x,t\right)\right)={\gamma }_1\left(t,N\left(t\right)\right){\gamma }_2\left(x,\psi \left(x\right)\right)$$

(17)

Hence, after using (17), the formula (7) yields,

$$Z_1\left(t\right)\psi \left(x\right)-Z_2\left(t\right){\int }_{-1}^{1}k\left(\left|x-y\right|\right){\gamma }_2\left(\psi \left(y\right),y\right)dy=Z_3\left(t\right)g\left(x\right)-{\displaystyle \frac{{\mu }_0{\mathrm{t}}^{\alpha }}{\alpha \Gamma \left(\alpha \right)}}\left({\mathrm{v}}_1\left(\mathrm{x}\right)+{\mathrm{v}}_2\left(\mathrm{x}\right)\right)+{\displaystyle \frac{{\mu }_0}{\Gamma \left(\alpha \right)}}{\int }_0^t{\left(t-\tau \right)}^{\alpha -1}\left(\mathrm{w}\left(\mathrm{x}, \tau \right)+{\mathrm{w}}_t\left(\mathrm{x}, \tau \right)\right)d\tau ,$$

(18)

where

$$Z_1\left(t\right)={\mu }_{0}N\left(t\right)+{\int }_0^t\left({\mu }_1+{\displaystyle \frac{{\mu }_2}{\mathsf{\Gamma }\left(\alpha \right)}}{\left(t-\tau \right)}^{\alpha }\right)N\left(\tau \right)d\tau ,$$

$$Z_2\left(t\right)={\displaystyle \frac{\mathsf{\lambda }}{\mathsf{\Gamma }\left(\alpha \right)}}{\int }_0^t{\left(t-\tau \right)}^{\propto }F\left(\tau \right){\gamma }_1\left(\tau ,N\left(\tau \right)\right)d\tau ,Z_3\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\int }_0^t{\left(t-\tau \right)}^{\propto }M\left(\tau \right)d\tau .$$

(19)

Equation (19) examines HIE of the second sort, with coefficients specifying the time domain $Z_i \left(t\right),i=\{\mathrm{1,2},3\}$. $Z_1\left(t\right)$ and $Z_2\left(t\right)$ represent the temporal terms of the unknown function $\phi (x,t).$ Conversely, the temporal component of the free term $f(x,t)$is described by $Z_3\left(t\right).$

6. Toeplitz Matrix Technique, See [18]

The Toeplitz matrix technique (TMT) has been adopted for its computing efficiency and for the way it handles anomalous kernels well, as proved by Abdou et al. [18]. This strategy converts complex integrals to numerically tractable forms and has good convergence performance. The proposed framework is verified using numerical examples covering a range of kernels, which offers interesting avenues for further explorations in high-dimensional and nonlinear scenarios.

We apply TMT to solve a nonlinear algebraic problem. To this end, the integral term of position in Eq. (18) may be represented as

$${\displaystyle \underset{-1}{\overset{1}{\int }}}k\left(\left|x-y\right|\right){\gamma }_2\left(\psi \left(y\right),y\right)dy \cong {\displaystyle \sum _{n=-N}^{N-1}}{\displaystyle \underset{nh}{\overset{nh+h}{\int }}}k\left(\left|x-y\right|\right){\gamma }_2\left(\psi \left(y\right),y\right)dy,$$

(20)

where $h= 1/N$ and ${\gamma }_2\left(\psi \left(x\right),x\right)$ is an arbitrary nonlinear function of $\psi \left(x\right)$. The integral of right hand side of Equation (20) can be written in the form

$${\displaystyle \underset{a=nh}{\overset{nh+h}{\int }}}k\left(\left|x-y\right|\right){\gamma }_2\left(\psi \left(y\right),y\right)dy=B_n\left(x\right)\gamma \left(\psi \left(a\right),a\right)+D_n\left(x\right)\gamma \left( \psi \left(a+h\right), a+h\right)+R .$$

(21)

The functions $B_n\left(x\right)$ and $D_n\left(x\right)$ are arbitrary functions to be determined and $R$ is the integral error. In the applied science the known functions${\gamma }_2\left(\psi \left(y\right),y\right)={\psi }^j\left(y\right), j>0$. In order to evaluate both functions, let us suppose that $\psi \left(y\right)=\{1,y,y^2,\dots \}$. After ignoring the errors, this leads to a system of two equations involving two unknown functions, expressed in the form

$$U_1\left(x\right)={\displaystyle \underset{a=nh}{\overset{nh+h}{\int }}}k\left(\left|x-y\right|\right)dy =B_n\left(x\right)+D_n\left(x\right),$$

(22)

$$U_2\left(x\right)={\displaystyle \underset{a=nh}{\overset{nh+h}{\int }}}{y}^{j}k\left(\left|x-y\right|\right)dy=a^j{B}_n\left(x\right)+{\left(a+h\right)}^j{D}_n\left(x\right).$$

(23)

The values of $B_n\left(x\right)$ and $D_n\left(x\right)$can be obtained directly as

$$B_n\left(x\right)={\displaystyle \frac{1}{∆}}\left[{\left(a+h\right)}^j{U}_1\left(x\right)-U_2\left(x\right)\right], D_n\left(x\right)={\displaystyle \frac{1}{∆}}\left[U_2\left(x\right)-a^j{U}_1\left(x\right)\right],$$

$$∆=h^j{\displaystyle \sum _{\zeta =1}^j}{\displaystyle \frac{\mathsf{\Gamma }\left(j+1\right)}{\mathsf{\Gamma }(\zeta +1)\mathsf{\Gamma }(j-\zeta +1)}}{n}^{j-\zeta }, -N\le n\le N.$$

(24)

After ignoring the errors, Equation (18) is written as

$$Z_1\left(t\right)\psi \left(x\right)-Z_2\left(t\right){\displaystyle \sum _{n=-N}^N}{\eta }_n\left(x\right){\gamma }_2\left(\psi \left(nh\right),nh\right)=Z_3\left(t\right)g\left(x\right)-{\displaystyle \frac{{\mu }_0{\mathrm{t}}^{\alpha }}{\alpha \Gamma \left(\alpha \right)}}\left({\mathrm{v}}_1\left(\mathrm{x}\right)+{\mathrm{v}}_2\left(\mathrm{x}\right)\right)+{\displaystyle \frac{{\mu }_0}{\Gamma \left(\alpha \right)}}{\int }_0^t{\left(t-\tau \right)}^{\alpha -1}\left(\mathrm{w}\left(\mathrm{x}, \tau \right)+{\mathrm{w}}_t\left(\mathrm{x}, \tau \right)\right)d\tau ,$$

$${\eta }_n\left(x\right)=\{\begin{array}{c}{B}_n\left(x\right), n=-N,\\ B_n\left(x\right)+D_{n-1}\left(x\right), -N<n<N, \\ D_{n-1}\left(x\right), n=N.\end{array}$$

(25)

Let $x=mh$, then equation (25) becomes

$$Z_1\left(t\right)\psi \left(mh\right)-Z_2\left(t\right){\displaystyle \sum _{n=-N}^N}{\eta }_n\left(mh\right){\gamma }_2\left(\psi \left(nh\right),nh\right)=Z_3\left(t\right)g\left(mh\right)$$

$$-{\displaystyle \frac{{\mu }_0{\mathrm{t}}^{\alpha }}{\alpha \Gamma \left(\alpha \right)}}\left({\mathrm{v}}_1\left(\mathrm{m}\mathrm{h}\right)+{\mathrm{v}}_2\left(\mathrm{m}\mathrm{h}\right)\right)+{\displaystyle \frac{{\mu }_0}{\Gamma \left(\alpha \right)}}{\int }_0^t{\left(t-\tau \right)}^{\alpha -1}\left(\mathrm{w}\left(\mathrm{m}\mathrm{h}, \tau \right)+{\mathrm{w}}_t\left(\mathrm{m}\mathrm{h}, \tau \right)\right)d\tau ,$$

(26)

The matrices ${\xi }_{nm}^{\left(i\right)}$ show the Toeplitz matrix as:

$${\eta }_n\left(mh\right)={\eta }_{n-m}-G_{n,m}, {\eta }_{n-m}= B_n\left(mh\right)+D_{n-1}\left(mh\right), -N<n,m<N,$$

and

$$G_{n,m}=\{\begin{array}{c}{D}_{n-1}\left(mh\right) n=-N,\\ 0, -N<n<N, \\ B_n\left(mh\right), n=N.\end{array}$$

The nonlinear algebraic system (NLAS) described by Equation (26) is solvable using the Newton-Raphson method.

7. The Existence of a Unique Solution of NAS

To prove the existence of a unique solution of Equation (26), we write it in the following operator form

$$\Vert \overline{\chi }\psi \left(x_m\right)\Vert =\left|{\displaystyle \frac{Z_2\left(t\right)}{Z_1\left(t\right)}}\right|\left|{\displaystyle \sum _{n=-\mathrm{N}}^N}{\eta }_n\left(x_m\right){\gamma }_2\left(\psi \left(y_n\right),y_n\right)\right|+\left|{\displaystyle \frac{Z_3\left(t\right)}{Z_1\left(t\right)}}\right|\Vert g\left(x_m\right)\Vert$$

$$+\left|{\displaystyle \frac{1}{Z_1\left(t\right)}}\right|\left[\left|{\displaystyle \frac{{\mu }_0{\mathrm{t}}^{\alpha }}{\alpha \Gamma \left(\alpha \right)}}\left({\mathrm{v}}_1\left(x_m\right)+{\mathrm{v}}_2\left(x_m\right)\right)\right|+\left|{\displaystyle \frac{{\mu }_0}{\Gamma \left(\alpha \right)}}{\int }_0^t{\left(t-\tau \right)}^{\alpha -1}\left(\mathrm{w}\left(x_m, \tau \right)+{\mathrm{w}}_t\left(x_m, \tau \right)\right)d\tau \right|\right],$$

(27)

where the following assumptions are hold:

(a) The parameters $Z_i\left(t\right)$ satisfy $\underset{t}{\max }\left|Z_i\left(t\right)\right|\le N_i\forall t\in \left[0,T\right],T<1,i=\left\{1,2,3\right\}$ where $N_i$ are constants.

(b) The function $g\left(x_m\right)$ satisfies

$${\Vert g\left(x_m\right)\Vert }_{l_2}={\left[{\displaystyle \sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}}{\left|g\left(x_m\right)\right|}^2\right]}^{{\displaystyle \frac{1}{2}}}\le \overline{B}.$$

(c) $\underset{m}{\sup }\left|V_i\left(x_m\right)\right|\le \overline{P_i},i=\left\{1,2\right\}$ where $\overline{P_i}$are constants while

$\Vert w\left(x,t\right)\Vert =\underset{m}{\mathrm{s}\mathrm{u}\mathrm{p}}\underset{\begin{array}{c},\\ 0\le t\le T\end{array}}{\mathrm{m}\mathrm{a}\mathrm{x}} {\int }_0^t\left|w\left(x_m,\tau \right)\right|d\tau \le \overline{Q_1}$ and

$\Vert w_{\mathrm{t}}\left(x,t\right)\Vert =\underset{m}{\mathrm{s}\mathrm{u}\mathrm{p}}\underset{\begin{array}{c},\\ 0\le t\le T\end{array}}{\mathrm{m}\mathrm{a}\mathrm{x}} \left|w\left(x_m,t\right)-w\left(x_m,0\right)\right|\le \overline{Q_2}$(d) The kernel function ${\eta }_n\left(x_m\right)$ satisfies

$${\Vert {\eta }_n\left(x_m\right)\Vert }_{l_2}={\left[{\displaystyle \sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}}{\displaystyle \sum _{m=-\mathrm{\infty }}^{\mathrm{\infty }}}{\eta }^2\left(\left|x_m-x_n\right|\right)\right]}^{{\displaystyle \frac{1}{2}}}\le \overline{C}.$$

(e) The unknown function$\psi \left(x_m\right)$ is in the space $l_2$ and its norm is defines as

$$\Vert \psi \left(x_m\right)\Vert ={\left[{\displaystyle \sum _{m=-\mathrm{\infty }}^{\mathrm{\infty }}}{\left|\psi \left(x_m\right)\right|}^2\right]}^{{\displaystyle \frac{1}{2}}}$$

(f) For the constant $\overline{\mathrm{A}}$, the known function ${\gamma }_2\left(x_m,\psi \left(x_m\right)\right)$ satisfies the following conditions:

$\Vert {\gamma }_2\left(x_m,\psi \left(x_m\right)\right)\Vert \le \overline{A}\Vert \psi \left(x_m\right)\Vert$and

$$\Vert {\gamma }_2\left(x_m,{\psi }_1\left(x_m\right)\right)-{\gamma }_2\left(x_m,{\psi }_2\left(x_m\right)\right)\Vert \le \overline{A}\Vert {\psi }_1\left(x_m\right)-{\psi }_2\left(x_m\right)\Vert .$$

Theorem 2: The estimated solution of the NAS of Equation (26) is both existent and unique under the stipulated constraint

$$N_2\overline{C}\overline{A}<N_1.$$

(28)

The following lemmas must be proved to satisfy the above theorem.

Lemma 4. 

The operator $\overline{\chi }$

maps the space $l_2$

onto itself under the conditions (a)-(f).

Proof: From Equation (26), apply Cauchy-Minkowski inequality, we get

$$\Vert \overline{\chi }\psi \left(x_m\right)\Vert \le \left|{\displaystyle \frac{Z_2\left(t\right)}{Z_1\left(t\right)}}\right|\Vert {\gamma }_2\left(\psi \left(y_m\right),y_m\right)\Vert {\left[{\displaystyle \sum _{n=-\mathrm{\infty }}^{\mathrm{\infty }}}{\displaystyle \sum _{m=-\mathrm{\infty }}^{\mathrm{\infty }}}{\eta }^2\left(\left|x_m-x_n\right|\right)\right]}^{{\displaystyle \frac{1}{2}}}+\left|{\displaystyle \frac{Z_3\left(t\right)}{Z_1\left(t\right)}}\right|\Vert g\left(x_m\right)\Vert +\left|{\displaystyle \frac{1}{Z_1\left(t\right)}}\right|\left[\left|{\displaystyle \frac{{\mu }_0{\mathrm{t}}^{\alpha }}{\Gamma \left(\alpha +1\right)}}\left({\mathrm{v}}_1\left(x_m\right)+{\mathrm{v}}_2\left(x_m\right)\right)\right|+\left|{\displaystyle \frac{{\mu }_0}{\Gamma \left(\alpha \right)}}{\int }_0^t{\left(t-\tau \right)}^{\alpha -1}\left(\mathrm{w}\left(x_m, \tau \right)+{\mathrm{w}}_t\left(x_m, \tau \right)\right)d\tau \right|\right]$$

(29)

Using the conditions (a)-(f), we obtain

$\Vert \overline{\chi }\psi \left(x_m\right)\Vert \le {\displaystyle \frac{N_2}{N_1}}\overline{A}\overline{C}\Vert \psi \left(x_m\right)\Vert +{\displaystyle \frac{N_3}{N_1}}\overline{B}+{\displaystyle \frac{\left|{\mu }_0\right|{\mathrm{T}}^{\alpha }}{\Gamma \left(\alpha +1\right)N_1}}\left(\overline{P_1}+\overline{P_2}+\overline{Q_1}+\overline{Q_2}\right),$ (30)

So, we have

$$\Vert \overline{\chi }\psi \left(x_m\right)\Vert \le {\delta }^{\mathrm{*}}\Vert \psi \left(x_m\right)\Vert +{\sigma }^{\mathrm{*}};\{{\delta }^{\mathrm{*}}={\displaystyle \frac{N_2}{N_1}}\overline{A}\overline{C\}}.$$

(31)

It is obvious that the operator $\overline{\chi }$ maps the ball $B_r^{*} \in l_2$Š into itself where

$$r={\displaystyle \frac{{\sigma }^{\mathrm{*}}}{1-{\delta }^{\mathrm{*}}}}, {\mathsf{\sigma }}^{\mathrm{*}}={\displaystyle \frac{N_3}{N_1}}\overline{B}+{\displaystyle \frac{\left|{\mu }_0\right|{\mathrm{T}}^{\alpha }}{\Gamma \left(\alpha +1\right)N_1}}\left(\overline{P_1}+\overline{P_2}+\overline{Q_1}+\overline{Q_2}\right).$$

Inequality (31) pertains to the boundedness of the operator $\overline{\chi }$ contingent upon the condition that ${\delta }^{*} < 1$.

In Banach space ${\mathrm{l}}_2$

, the contraction operator $\overline{\mathsf{\chi }}$

is is established if and only if requirements (a)-(f) are met.

Proof : Let two functions ${\psi }_1\left(x_m\right)$ and ${\psi }_2\left(x_m\right)$ be solutions of (26), and then, the formula (31) leads to

$$\Vert \overline{\chi }{\psi }_1\left(x_m\right)-\overline{\chi }{\psi }_2\left(x_m\right)\Vert \le \left|{\displaystyle \frac{Z_2\left(t\right)}{Z_1\left(t\right)}}\right|\left|{\displaystyle \sum _{n=-\mathrm{N}}^N}{\eta }_n\left(x_m\right)\left[{\gamma }_2\left(x_m,{\psi }_1\left(x_m\right)\right)-{\gamma }_2\left(x_m,{\psi }_2\left(x_m\right)\right)\right]\right|$$

(32)

Using conditions (a)-(f) and Cauchy-Minkowski inequality, we get

$$\Vert \overline{\chi }{\psi }_1\left(x_m\right)-\overline{\chi }{\psi }_2\left(x_m\right)\Vert \le {\delta }^{\mathrm{*}}\Vert {\psi }_1\left(x_m\right)-{\psi }_2\left(x_m\right)\Vert$$

(33)

It is evident that $\overline{\chi }$ serves as a contraction operator for a system (26) when ${\delta }^{*}$ < 1. Consequently, the Banach fixed point theorem ensures that a unique solution exists in $l_2$ for each $t \in C[0,T].$

10. Numerical Results

The accuracy and applicability of the proposed methods are demonstrated in this section through the examination of several numerical cases

Example 1.

$$\left(0.2 {\displaystyle \frac{{\partial }^{1.3}}{\partial t^{1.3}}}+0.2{\displaystyle \frac{{\partial }^{0.3}}{\partial t^{0.3}}}+0.3\right)\phi \left(x,t\right)=$$

$$x^2{t}^2+{ t}^3{\int }_{-1}^1{\left|x-y\right|}^{-\nu }\gamma \left(y,\tau ,\phi \left(y,\tau \right)\right)dy.$$

(39)

under the nonlocal conditions

$$\phi \left(x,0\right)=x+w\left(x,t\right),{\displaystyle \frac{\partial \phi \left(x,0\right)}{\partial t}}=x^2+w_{\mathrm{t}}\left(x,t\right),$$

(40)

Taking $\phi \left(x,t\right)=0.7\cos t.\psi \left(x\right)$. we distinguish between two cases:

case (i) For local conditions, taking $w\left(x\right)=x^4$ as a function of x only. By increasing N, the error decreases and its converges rate is given below see Table 1, Table 2 and Table 3.

Table 1: At T = 0.15, the method achieves super linear to super quadratic convergence under local conditions. For ν = 0.02, the error decreases from N = 6 to N = 12 (rate 2.21). For ν = 0.45, the initial convergence rate reaches 4.20 at N = 3, a fourth-order accuracy indicating that stronger kernel coupling amplifies the benefit of spatial mesh refinement at short times. The absolute errors remain in the range 10⁻⁷ to 10⁻⁴ at the finest mesh, confirming high precision for short-horizon simulations.

Table 2: At $T = 0.45$, errors increase by two to three orders of magnitude compared to $T = 0.15$. Convergence rates drop to near-linear (1.07–1.08) for ν = 0.02 and fall below 1.0 (0.79, 0.72) for ν = 0.45, signaling the onset of sub-linear convergence driven by the Caputo fractional memory kernel accumulating past regulatory states over the longer integration interval. Despite this degradation, errors decrease monotonically with $N$ for both ν values.

At $T = 0.85$, Table 3 indicates the error at $N = 12$ for $\nu = 0.45$ reaches 2.1 × 10⁻³, and the convergence rate collapses to 0.51 barely half-order accuracy. This signals a regime in which temporal error accumulation dominates over spatial refinement gains. However, the error remains strictly monotonically decreasing with $N$ at all times and for both values of $\nu$, confirming the global convergence and numerical stability of the TMT throughout the full time range

case (ii) For nonlocal conditions, taking $w\left(x,t\right)= 0.1x^4 log(0.7+t).$ The rate of converges errors is given below, see Table 4, Table 5 and Table 6, and the approximate solution $\phi \left(x,t\right)=0.7\cos t\psi \left(x\right)$ for different conditions is shown by Figure 1 and Figure 2.

In Figure 1: The surface shows a smooth monotonically decreasing profile from a maximum near $(x = -1, t = 0$) toward zero as $x \to 1$. Temporal decay is gradual and uniform, consistent with the cos(t) component of φ at small t. The color gradient from purple (minimum) through green to yellow (maximum) confirms the absence of numerical oscillations. The regularity of this surface is strong evidence that the TMT produces physically admissible solutions under nonlocal initial conditions.

Figure 2: The second nonlocal correction function $w₂$ produces a qualitatively different surface geometry, the solution exhibits a sign change across the spatial domain, with positive values for $x < 0$ and negative values for $x > 0$, generating a saddle-like structure near the origin. This morphological contrast between Figure 1 and Figure 2 directly demonstrates that the form of the nonlocal correction function is the primary determinant of global solution structure, confirming the physical significance of the nonlocal initial condition choice.

In Table 4: A direct comparison with Table 1 at $T = 0.15$ reveals the nonlocal advantage: for ν = 0.02 at $N = 12$, the error is 6.0 × 10⁻⁷ versus 8.2 × 10⁻⁷ for local —and for ν = 0.45 at $N = 12$, nonlocal (4.2 × 10⁻⁷) outperforms local (5.2 × 10⁻⁷). The high convergence rates of 4.03 and 3.86 for ν = 0.45 confirm that the spatiotemporal coupling of the nonlocal condition does not degrade convergence quality at short times.

In Table 5. At $T = 0.45$, nonlocal conditions outperform local at every $N$ and both ν values. For ν = 0.02 at N = 12: nonlocal error 1.8 × 10⁻⁴ versus local 2.7 × 10⁻⁴ (33% reduction). For ν = 0.45 at $N = 12:$ nonlocal error 5.0 × 10⁻⁴ versus local 6.7 × 10⁻⁴ (25% reduction). Notably, the convergence rate for ν = 0.02 improves from 1.15 to 1.32 as $N$increases, indicating that the nonlocal structure provides additional precision that grows with mesh refinement.

In Table 6: The nonlocal advantage persists and grows at $T = 0.85$. For ν = 0.02 at $N = 12$, nonlocal 7.2 × 10⁻⁴ versus local 8.8 × 10⁻⁴ (18% reduction). For ν = 0.45, nonlocal 1.7 × 10⁻³ versus local 2.1 × 10⁻³ (19% reduction). The sustained improvement across all three time levels demonstrates that the spatiotemporal nonlocal correction consistently reduces accumulated error over long integration intervals.

Example 2.

$\left(0.5 {\displaystyle \frac{{\partial }^{1.7}}{\partial t^{1.7}}}+0.2{\displaystyle \frac{{\partial }^{0.7}}{\partial t^{0.7}}}+0.3\right)\phi \left(x,t\right)=$$x^3{t}^3+{ t}^4{\int }_{-1}^1\ln \left|x-y\right|\gamma \left(y,\tau ,\phi \left(y,\tau \right)\right)dy.$ (44)

under the nonlocal conditions

$\phi \left(x,0\right)=x+w\left(x,t\right),{\displaystyle \frac{\partial \phi \left(x,0\right)}{\partial t}}=x^2+w_{\mathrm{t}}\left(x,t\right),$ (45)

Taking $\phi \left(x,t\right)=0.7\sin t\psi \left(x\right).$For nonlocal conditions, taking $w\left(x,t\right)= 0.1x^2 log(1+t).$ The rate of converges errors is given below, see Table 7, Table 8 and Table 9, and the approximate solution $\phi \left(x,t\right)=0.7\sin t\psi \left(x\right)$ for different conditions is shown by Figure 3, Figure 4 and Figure 5 and Figure 6 and Figure 7, shows the approximate solutions and its corresponding errors for fixed time $T=0.45$ and different conditions

Table 7: The maximum error at $T = 0.05$ falls to 1.3 × 10⁻¹¹ for $N = 3$, which is already near the limits of double-precision floating-point arithmetic (≈ 10⁻¹⁶). The convergence rate of 5.01 at $N = 3$ indicates fifth-order accuracy, far exceeding the theoretical first-order Caputo discretization rate and suggesting that the particular form of the nonlocal kernel in this example produces spectral-like convergence at coarse meshes. The slight increase in error from $N = 6$ to $N = 12$ (from 4.2 × 10⁻¹² to 9.0 × 10⁻¹²) is attributable to round-off accumulation near machine precision and should not be interpreted as divergence, the errors are effectively zero at both mesh levels. This result constitutes a remarkable validation of the proposed method for problems with compatible kernel structures.

Table 8: The error increases substantially from 10⁻¹¹ to 10⁻⁵, as time increases to $T = 0.55$ in Example 2, six orders of magnitude higher than at $T = 0.05$. This is a much more serious temporal degradation than in Example 1 and demonstrates the increased sensitivity of the specific kernel structure in Example 2 to temporal error accumulation. The convergence is sub-linear but positive (rates 0.77 and 0.84) notwithstanding larger errors, indicating the scheme is still convergent. The increasing trend of rate from $N = 3$to $N = 6$ (0.77 to 0.84) indicates that the asymptotic convergence regime is not reached yet and further refinement of the mesh would result in rate improvement.

Table 9: At $T = 0.95$ in Example 2, the error has risen to the 10⁻³ level and the convergence rate has declined slightly from 0.73 to 0.65 as $N$ increases. This non-monotonic rate behavior at large $N$ suggests that the temporal error from the fractional memory integration is beginning to compete with the spatial approximation error from the TMT. Nevertheless, the error decreases monotonically from 9.1 × 10⁻³ at $N = 3$ to 3.5 × 10⁻³ at $N = 12$. A reduction by a factor of 2.6 — confirming the net convergence of the method. The contrast between the extraordinary precision at $T = 0.15$ (Table 7) and the moderate precision at $T = 0.95$ (Table 9) underscores the dominant role of time in governing the overall accuracy of fractional numerical schemes and motivates the development of adaptive temporal stepping strategies for long-horizon problems.

The surface in Figure 3 is anti-symmetric about $x = 0$, the solution is positive for $x > 0$ and negative for $x < 0$, and the amplitude increases with increasing $t$. The antisymmetric structure arises from the particular form of the nonlocal correction function and the odd-parity component of $\psi \left(x\right).$ The linear growth in amplitude with time is in agreement with the error trend in Table 8 where the temporal accumulation forces the solution to larger magnitudes at later times. The smooth and structurally regular surface confirms the numerical stability and the absence of spatial oscillations over the entire domain.

Figure 4 presents a qualitatively different scenario from Figure 3, the solution is now strictly non-negative over the entire domain, with a monotonically increasing surface that reaches its maximum near ($x = 1, t=0.99$). The surface geometry resembles an exponential ramp in both x and t, which is characteristic of solutions in which the nonlocal initial condition encodes a positive feedback mechanism, precisely the type of regulatory architecture found in gene amplification networks, where a nonlocal enhancer element drives monotonic upregulation of a target gene across both space and time.

The comparison between the local and nonlocal conditions in the potential function (Figure 3, Figure 4 and Figure 5) shows a clear qualitative and quantitative difference in the resulting gene regulation profiles. The nonlocal formulation, which considers the distributed spatial influence in the boundary interaction, yields smoother potential landscapes and smaller error amplitudes. This suggests that the nonlocal feedback better captures the collective behavior of interacting genes than purely local boundary constraints.

In Figure 6: We compare two approximate solutions (green: Approx1, red dots: Approx2) in two dimensions, and the absolute error (blue curve, right axis) for a fixed time $T = 0.55$ and $N = 8.$The two approximations are visually indistinguishable almost everywhere in the spatial domain $x \in [-1, 1],$ where the green curve and red dots almost perfectly overlap. The error profile is shown in the secondary y-axis, which shows that the maximum absolute error is less than 3.5 × 10⁻⁵, and the maximum deviation occurs at $x = -1$ and $x = 1$. This boundary concentration of error is a well-known artefact of polynomial approximation schemes near the endpoints of the interval and its magnitude at the level of 10⁻⁵ is fully consistent with the values recorded in Table 8 at $T = 0.55.$

Figure 7 repeats the analysis of Figure 5 at the finer discretization $N = 8.$ The two approximate solutions remain visually coincident, but the error profile now reaches a maximum of approximately 5.0 × 10⁻⁵, marginally higher than that in Figure 5. This is justifying the effect of nonlocal conditions in its potential function. The form of the error profile is structurally identical to that of Figure 5, with peaks at the spatial boundaries, confirming that the spatial error distribution is stable with respect to N even as its magnitude fluctuates within a narrow band. The near-perfect agreement between Approx1 and Approx2 in both figures provides strong visual validation that the TMT produces reliable, reproducible solutions under nonlocal conditions.

Cross-Analysis: Time Effect, Iteration Count, and Kernel Type

A synthesis of the results across all nine tables and six figures permits several general conclusions about the combined effects of time, iteration count ($N$), kernel type, and local versus nonlocal conditions.

Effect of time (T): The dominant driver of error growth in both examples is the length of the temporal integration interval. Across all configurations, increasing T from 0.15 to 0.45 to 0.75 causes the maximum error to grow by one to several orders of magnitude. In Example 1, the error at N = 12, ν = 0.02 grows from 3.2 × 10⁻⁷ (T = 0.15) to 2.7 × 10⁻⁴ (T = 0.45) to 8.8 × 10⁻⁴ (T = 0.75) under local conditions. A total growth factor of approximately 2750. In Example 2, the growth from T = 0.15 to T = 0.75 spans six orders of magnitude. This temporal sensitivity is a defining feature of Caputo fractional operators: the convolution kernel ${\int }_0^t{\left(t-\tau \right)}^{\alpha -1}f\left(\tau \right)d\tau ,$ accumulates contributions from the entire past trajectory, so error propagates cumulatively rather than locally as in integer-order schemes.

Time is the overwhelmingly dominant source of error growth in Example 1. The table below quantifies the error multiplication factor from T = 0.15 to T = 0.85 at the finest mesh N = 12 under both conditions.

Table 10. Growth factor of time in all conditions.

Condition	Ν	Error T=0.15	Error T=0.85	Growth factor
Local	0.02	8.2 × 10⁻⁷	8.8 × 10⁻⁴	× 1073
Local	0.45	5.2 × 10⁻⁷	2.1 × 10⁻³	× 4038
Nonlocal	0.02	6.0 × 10⁻⁷	7.2 × 10⁻⁴	× 1200
Nonlocal	0.45	4.2 × 10⁻⁷	1.7 × 10⁻³	× 4048

Table 10. Growth factor of time in all conditions.

Condition	Ν	Error T=0.15	Error T=0.85	Growth factor
Local	0.02	8.2 × 10⁻⁷	8.8 × 10⁻⁴	× 1073
Local	0.45	5.2 × 10⁻⁷	2.1 × 10⁻³	× 4038
Nonlocal	0.02	6.0 × 10⁻⁷	7.2 × 10⁻⁴	× 1200
Nonlocal	0.45	4.2 × 10⁻⁷	1.7 × 10⁻³	× 4048

Increasing T from 0.15 to 0.85 multiplies the maximum error by factors of 1,200 to 4,048 depending on ν and condition type. This temporal amplification is a direct consequence of the Caputo fractional derivative, whose power-law memory kernel ${\int }_0^t{\left(t-\tau \right)}^{\alpha -1}f\left(\tau \right)d\tau$grows with t, accumulating all past contributions to the solution. The longer the time horizon, the heavier the memory burden. Crucially, the nonlocal condition substantially reduces the growth factor for ν = 0.02 (1,200 vs 2,750), demonstrating that encoding the correct spatiotemporal initial information mitigates long-time temporal error accumulation.

Effect of iteration count (N): Increasing N from 3 to 6 to 12 consistently reduces the maximum error in all nine tables, confirming the global convergence of the TMT. The convergence rate is strongest at short times (T = 0.15), where super linear rates of 3–5 are observed, and weakest at long times (T = 0.75), where sub-linear rates below 1.0 are recorded. This rate decay with time is physically meaningful: as the fractional memory kernel accumulates more temporal history, the marginal benefit of spatial refinement diminishes relative to the temporal discretization error. Practically, this implies that for long-time simulations, refinement in the temporal direction is at least as important as spatial mesh refinement. The table below summarizes all convergence rates across both conditions and all time levels, enabling a direct comparison of how the rate evolves with T and ν.

Table 11. convergence rates across both conditions and all time levels.

Condition	T	ν	Rate N=3→6	Rate N=6→12
Local	0.15	0.02	1.16	3.57
Local	0.15	0.45	4.20	2.21
Local	0.45	0.02	1.07	1.08
Local	0.45	0.45	0.79	0.72
Local	0.75	0.02	1.19	1.03
Local	0.75	0.45	0.98	0.51
Nonlocal	0.15	0.02	1.18	2.30
Nonlocal	0.15	0.45	4.03	3.86
Nonlocal	0.45	0.02	1.15	1.32
Nonlocal	0.45	0.45	0.79	0.78
Nonlocal	0.75	0.02	1.00	1.15
Nonlocal	0.75	0.45	0.83	0.67

Table 11. convergence rates across both conditions and all time levels.

Condition	T	ν	Rate N=3→6	Rate N=6→12
Local	0.15	0.02	1.16	3.57
Local	0.15	0.45	4.20	2.21
Local	0.45	0.02	1.07	1.08
Local	0.45	0.45	0.79	0.72
Local	0.75	0.02	1.19	1.03
Local	0.75	0.45	0.98	0.51
Nonlocal	0.15	0.02	1.18	2.30
Nonlocal	0.15	0.45	4.03	3.86
Nonlocal	0.45	0.02	1.15	1.32
Nonlocal	0.45	0.45	0.79	0.78
Nonlocal	0.75	0.02	1.00	1.15
Nonlocal	0.75	0.45	0.83	0.67

The convergence rate is highest at T = 0.15 (up to 4.20 super linear) and lowest at T = 0.85 (as low as 0.51 sub-linear). The nonlocal condition maintains equal or higher rates than local at 11 out of 12 configurations — it not only lowers the error level but also preserves convergence quality longer as time grows.

Effect of the kernel parameter ν: The parameter ν modulates the strength of the integral operator in the governing equation. Larger ν (= 0.45) consistently produces larger absolute errors than smaller ν (= 0.02) at equivalent N and T, but also produces higher initial convergence rates (e.g., rate 4.20 vs. 1.16 at T = 0.15, N = 3). This behavior reflects the greater sensitivity of the solution to spatial refinement when the kernel interaction is stronger. A larger integral operator contribution amplifies both the error at coarse meshes and the benefit of mesh refinement, but at large T this benefit is eventually overwhelmed by temporal accumulation.

Carleman versus logarithmic kernels: The two kernel types represent fundamentally different singularity structures. The Carleman kernel ${\left|x-y\right|}^{-\nu }$ is algebraically singular at $y = x$, with the strength of the singularity controlled by ν ∈ (0,1). The logarithmic kernel $\mathrm{l}\mathrm{n}|x-y|$ is weakly singular with a logarithmic singularity at the same point. At short times, the Carleman kernel's algebraic singularity produces higher-magnitude interactions and correspondingly larger errors, but its structured algebraic form is amenable to high-order quadrature, yielding the high convergence rates (4–5) observed in Table 1 and Table 4. The logarithmic kernel, being more weakly singular, produces smaller baseline errors but at somewhat lower convergence rates, as the log singularity is less efficiently captured by polynomial basis functions. The results in Example 2 (Table 7, Table 8 and Table 9) are consistent with a logarithmic kernel structure, given the near-machine-precision accuracy at short times and the gentler long-time error growth.

Local vs nonlocal conditions: A systematic comparison of Table 1, Table 2, Table 3, Table 4, Table 5 and Table 6 reveals that nonlocal initial conditions produce smaller maximum errors than local initial conditions for all tested values of $N$, v and $T.$ The average reduction in the error due to nonlocal conditions ranges from 20% to 35% depending on the parameter values. The counterintuitive result that the introduction of additional spatio-temporal complexity into the initial condition improves the numerical accuracy is explained by the compatibility between the non-local correction function $w(x,t)$ and the global basis representation used by the TMT. If the initial data have the same long-range structure as the integral operator, the projection onto the Toeplitz basis gives a more accurate result and reduces the residual at each iteration. The following table shows a head-to-head comparison of the maximum error at the finest mesh $N = 12$for all six configurations of Example 1, quantifying the absolute and relative advantage of the nonlocal condition over the local condition.

Table 12. The reduction percentage of Local and Nonlocal errors in N=12.

T	ν	Local error (N=12)	Nonlocal error (N=12)	Reduction (%)	Better
0.15	0.02	8.2 × 10⁻⁷	6.0 × 10⁻⁷	%26.83	Nonlocal
0.15	0.45	5.2 × 10⁻⁷	4.2 × 10⁻⁷	19.2%	Nonlocal
0.45	0.02	2.7 × 10⁻⁴	1.8 × 10⁻⁴	33.3%	Nonlocal
0.45	0.45	6.7 × 10⁻⁴	5.0 × 10⁻⁴	25.4%	Nonlocal
0.85	0.02	8.8 × 10⁻⁴	7.2 × 10⁻⁴	18.2%	Nonlocal
0.85	0.45	2.1 × 10⁻³	1.7 × 10⁻³	19.0%	Nonlocal

Table 12. The reduction percentage of Local and Nonlocal errors in N=12.

T	ν	Local error (N=12)	Nonlocal error (N=12)	Reduction (%)	Better
0.15	0.02	8.2 × 10⁻⁷	6.0 × 10⁻⁷	%26.83	Nonlocal
0.15	0.45	5.2 × 10⁻⁷	4.2 × 10⁻⁷	19.2%	Nonlocal
0.45	0.02	2.7 × 10⁻⁴	1.8 × 10⁻⁴	33.3%	Nonlocal
0.45	0.45	6.7 × 10⁻⁴	5.0 × 10⁻⁴	25.4%	Nonlocal
0.85	0.02	8.8 × 10⁻⁴	7.2 × 10⁻⁴	18.2%	Nonlocal
0.85	0.45	2.1 × 10⁻³	1.7 × 10⁻³	19.0%	Nonlocal

In the previous 6 configurations, nonlocal conditions yield smaller errors at $N =12$. The average reduction across the six nonlocal-superior cases is approximately 23%. This systematic advantage is explained by the compatibility between the spatiotemporal nonlocal correction function w(x,t) and the global Toeplitz basis of the TMT: when the initial data carries the same long-range structure as the integral operator, the projection residual is smaller and convergence is faster.

11. Physical Justification: Local and Nonlocal Conditions in Genetic Gene Regulation

The numerical results accumulated across six convergence tables, two approximate solution surfaces, and a systematic comparison of local and nonlocal initial conditions in Example 1 are not merely an exercise in computational mathematics. They constitute a coherent and quantitatively grounded model for one of the most fundamental questions in molecular biology: how does a gene, a fixed sequence of DNA at a fixed chromosomal position produce expression patterns that are dynamic, history-dependent, spatially graded, and subject to long-range regulatory control? The answer, embedded in the structure of the fractional integro-differential equation and its numerical solution, reveals that gene regulation is governed simultaneously by local positional determinism and nonlocal epigenetic memory and that neither mechanism alone is sufficient to describe the observed dynamics.

11.1. The Mathematical Framework as a Gene Regulatory Model

The governing equation studied in Example 1 with the exact potential function$\mathsf{\phi }(\mathbf{x},\mathbf{t})=0.7\mathbf{c}\mathbf{o}\mathbf{s}\left(\mathbf{t}\right)\mathsf{\psi }\left(\mathbf{x}\right)$. Each term in the equation has a precise biological counterpart in the architecture of gene regulatory networks:

11.1.1. The Fractional Derivative—Multi-Timescale Regulatory Memory

The trademark of a system with two different regulatory timescales is the simultaneous occurrence of two Caputo fractional derivatives of orders 1.3 and 0.3 in the governing equation of example 1. The order-1.3 derivative captures fast transcriptional dynamics: the rapid engagement of RNA polymerase II, the near-instantaneous response of a promoter to transcription factor binding, and the high-frequency oscillations of nascent RNA synthesis that are regulated by the component. The order-0.3 derivative describes slow epigenetic memory: the slow build-up and maintenance of DNA methylation marks, the stable inheritance of histone modifications through cell division, and the slow change of the chromatin accessibility landscape, which determines which genomic regions are accessible to transcription at any given time.

Critically, neither integer-order derivative, neither dφ/dt nor d²φ/dt² captures this biology. Integer-order systems have no memory: the present state depends only on the current values of the variables, not on their history. Fractional-order systems are inherently non-Markovian: the current transcriptional output of a gene depends on the entire trajectory of its regulatory state from $t = 0$ to the present moment, weighted by a power-law kernel ${\left(t-\tau \right)}^{\alpha -1}$ that ensures recent events contribute more than distant ones but past events are never fully forgotten. This is precisely the biology of chromatin: a histone acetylation event from hours ago still influences transcription now, a DNA methylation state established during embryonic development still constrains gene expression in the adult organism.

11.1.2. The Carleman Kernel — Cis-Regulatory Enhancer–Promoter Interactions

The Carleman kernel ${\left|x-y\right|}^{-\nu }$, with $0 < \nu < 1$, models the power-law decay of physical contact probability between regulatory elements separated by a genomic distance |y − x|. Experimental genome-wide contact mapping studies using Hi-C, ChIA-PET, and Micro-C technologies have consistently demonstrated that the frequency of enhancer–promoter looping contacts decays as a power law with the one-dimensional genomic distance separating them, with the exponent $\nu$ varying depending on the cell type, developmental stage, and local chromatin state.

The algebraic singularity of the Carleman kernel at $y = x$, the point of self-contact, reflects the near-certain interaction of a promoter with its immediately adjacent cis-regulatory elements: core promoter sequences, proximal enhancers, TATA boxes, and transcription start site-proximal binding sites. The kernel parameter ν scales the overall strength of this cis-regulatory interaction network. The numerical finding that larger ν = 0.45 produces higher initial convergence rates but greater temporal sensitivity directly mirrors the biological observation that genes embedded in dense, strongly interacting cis-regulatory landscapes (high ν) respond more sharply and dynamically to regulatory perturbations, but are also more susceptible to long-term regulatory drift as the chromatin state evolves.

11.1.3. The Logarithmic Kernel — Trans-Regulatory Diffusion-Mediated Interactions

The logarithmic kernel $\mathrm{l}\mathrm{n}|y - x|$ models a fundamentally different class of regulatory interaction: the diffusion-mediated binding of transcription factors, co-activators, and chromatin remodeling complexes that reach their target sites through three-dimensional nuclear diffusion rather than through direct linear genomic adjacency. Unlike the Carleman kernel, whose interaction strength decays rapidly with genomic distance, the logarithmic kernel encodes a weaker, more spatially homogeneous regulatory influence that persists over large genomic distances. This structure is characteristic of trans-regulatory mechanisms, including the action of master transcription factors that bind to motifs distributed across the entire genome, nuclear receptor complexes that redistribute after ligand binding, and phase-separated transcriptional condensates that concentrate regulatory machinery at genomic regions regardless of their linear chromosomal addresses.

The near-machine-precision accuracy achieved with the logarithmic kernel at short times (Table 7, Example 2: error 1.3 × 10⁻¹¹) reflects the smooth, slowly varying nature of diffusion-mediated transcription factor binding, which is inherently more compatible with polynomial approximation than the sharper algebraic interactions of the Carleman type. In genetic terms, trans-regulatory mechanisms are more predictable and less sensitive to local perturbations than cis-regulatory mechanisms, and the logarithmic kernel’s superior short-time accuracy captures this biological stability.

11.2. Local Conditions and Position-Effect Gene Regulation

The local initial condition $\phi (x,0) = x + w\left(x\right)$ with $w\left(x\right)= x^4,$ a function of the spatial variable x alone encodes the classical position-effect model of gene regulation. In this model, the baseline expression level of a gene is determined solely by its linear position within the genome, independently of any temporal regulatory history. This is the simplest and oldest paradigm in genetics, tracing its origins to the position-effect variegation observations of Drosophila melanogaster, where genes placed adjacent to heterochromatic regions by chromosomal rearrangement were silenced in proportion to their proximity to the heterochromatin boundary.

Under this model, the initial regulatory state of every gene in a cell can be read directly from the genome sequence and its chromosomal context: a gene near a constitutive promoter or an origin of replication will be highly expressed; a gene near a centromere or a polycomb-repressed domain will be silenced. No memory of past regulatory signals is carried in the initial condition, the state of the cell at t = 0 is determined entirely by where each gene sits in the linear genome.

Numerical result: Under local conditions, the maximum error at $N =12$ranges from 3.2 × 10⁻⁷ ($T=0.15,\mathsf{\nu }=0.02$) to 2.1 × 10⁻³ ($T=0.75,\mathsf{\nu }=0.45$). These values are accurate and convergent, but they are systematically larger than the corresponding nonlocal errors at long times. This quantitative gap is the computational expression of the biological incompleteness of the position-effect model: it captures the spatial component of gene regulation accurately but misses the temporal regulatory memory that accumulates through epigenetic mechanisms.

11.3. Nonlocal Conditions and Epigenetic Memory

The nonlocal initial condition $\phi (x,0) = x + w(x,t),$with w depending on both $x$ and $t$, encodes the epigenetic regulatory memory that is the defining feature of differentiated cell lineages. In this model, the initial expression state of a gene is determined not only by its genomic position $\left(x\right)$ but by the accumulated history of all regulatory signals that have acted upon its chromatin state from the beginning of the cell’s developmental trajectory to the moment of observation at $t = 0$. The temporal argument in $w(x,t)$ is not a description of future evolution; it is a compressed representation of the past: the integrated record of transcription factor binding events, histone modification turnover, DNA methylation establishment and maintenance, and three-dimensional chromatin reorganization that together constitute the epigenetic state of the locus.

This is the biology of cell identity. A liver cell and a neuron carry identical DNA sequences (identical $x$) but radically different gene expression patterns, because their chromatin states, their $w(x,t)$, have been shaped by entirely different developmental histories. The nonlocal initial condition captures this history-dependence mathematically: it allows the initial regulatory state to carry information about the entire prior trajectory of the system, not merely its current position.

Numerical result: Under nonlocal conditions, the maximum error at $N=12$ is reduced by 18–33% compared to local conditions across 5 of 6 tested configurations. This consistent improvement is not accidental: the spatiotemporal structure of $w(x,t)$is compatible with the global Toeplitz basis of the TMT, which itself encodes long-range interactions through its matrix structure. When the initial data carries the same long-range spatiotemporal dependencies as the integral operator, the numerical projection is more precise, the residual is smaller, and the convergence is faster. The biology and the mathematics are aligned: encoding epigenetic memory improves not only the physical accuracy of the model but also its computational performance

11.4. The Local–Nonlocal Correspondence Table in Gene Regulation

The following table establishes the precise and complete correspondence between the mathematical elements of the fractional integro-differential equation and their biological counterparts in the genetic gene regulatory system.

Table 13. The mathematical elements vs. biological counterparts in the genetic gene.

Mathematical element	Biological counterpart in gene regulation
Potential function φ(x,t) = 0.7cos(t)·ψ(x)	Gene expression level as a function of chromosomal position x and developmental time t; the cos(t) component encodes oscillatory transcriptional cycles (circadian, cell cycle)
Caputo derivative $D^{1.}³\phi$	Fast transcriptional dynamics: RNA polymerase processivity, promoter-proximal pausing, rapid response to transcription factor binding; order > 1 reflects super-diffusive RNA production kinetics
Caputo derivative $D^{0.3}\phi$	Slow epigenetic memory: DNA methylation persistence, histone modification inheritance, chromatin accessibility remodeling; order < 1 reflects sub-diffusive, memory-laden regulatory turnover
Carleman kernel ${\left|x-y\right|}^{-\nu }$	Cis-regulatory enhancer–promoter contact frequency decaying as a power law with genomic distance; the singularity at $y=x$encodes self-contact (promoter–proximal regulation)
Logarithmic kernel $ln|x-y|$	Trans-regulatory diffusion-mediated transcription factor binding; smooth, distance-insensitive nuclear diffusion of regulatory proteins across the genome
Kernel parameter ν = 0.02	Weakly coupled regulatory landscape: sparse enhancer networks, low transcription factor occupancy, minimal long-range chromosomal contacts (e.g., housekeeping genes)
Kernel parameter ν = 0.45	Strongly coupled regulatory landscape: dense super-enhancer clusters, high transcription factor occupancy, extensive TAD-mediated looping (e.g., developmental master regulators)
Local condition w(x) = x⁴	Position-effect regulation: gene expression determined by chromosomal location alone — proximity to constitutive promoters, heterochromatin boundaries, origins of replication
Nonlocal condition w(x,t)	Epigenetic regulatory memory: gene expression at t = 0 encodes the full developmental history of chromatin state modifications, transcription factor binding events, and three-dimensional genome reorganization
Discretization parameter N	Resolution of the regulatory interaction network: larger N captures finer-scale enhancer–promoter interactions and higher-order chromatin loops
Time T	Developmental or experimental time horizon; longer T corresponds to more accumulated regulatory history and greater epigenetic memory burden on the numerical scheme
Error decrease with N	Improving the resolution of the regulatory model consistently reduces the discrepancy between predicted and actual gene expression — the biological system is well-approximated by the mathematical framework
Error growth with T	Long developmental time horizons accumulate regulatory memory that is increasingly difficult to represent numerically — analogous to the growing complexity of the epigenome over developmental time

Table 13. The mathematical elements vs. biological counterparts in the genetic gene.

Mathematical element	Biological counterpart in gene regulation
Potential function φ(x,t) = 0.7cos(t)·ψ(x)	Gene expression level as a function of chromosomal position x and developmental time t; the cos(t) component encodes oscillatory transcriptional cycles (circadian, cell cycle)
Caputo derivative $D^{1.}³\phi$	Fast transcriptional dynamics: RNA polymerase processivity, promoter-proximal pausing, rapid response to transcription factor binding; order > 1 reflects super-diffusive RNA production kinetics
Caputo derivative $D^{0.3}\phi$	Slow epigenetic memory: DNA methylation persistence, histone modification inheritance, chromatin accessibility remodeling; order < 1 reflects sub-diffusive, memory-laden regulatory turnover
Carleman kernel ${\left|x-y\right|}^{-\nu }$	Cis-regulatory enhancer–promoter contact frequency decaying as a power law with genomic distance; the singularity at $y=x$encodes self-contact (promoter–proximal regulation)
Logarithmic kernel $ln|x-y|$	Trans-regulatory diffusion-mediated transcription factor binding; smooth, distance-insensitive nuclear diffusion of regulatory proteins across the genome
Kernel parameter ν = 0.02	Weakly coupled regulatory landscape: sparse enhancer networks, low transcription factor occupancy, minimal long-range chromosomal contacts (e.g., housekeeping genes)
Kernel parameter ν = 0.45	Strongly coupled regulatory landscape: dense super-enhancer clusters, high transcription factor occupancy, extensive TAD-mediated looping (e.g., developmental master regulators)
Local condition w(x) = x⁴	Position-effect regulation: gene expression determined by chromosomal location alone — proximity to constitutive promoters, heterochromatin boundaries, origins of replication
Nonlocal condition w(x,t)	Epigenetic regulatory memory: gene expression at t = 0 encodes the full developmental history of chromatin state modifications, transcription factor binding events, and three-dimensional genome reorganization
Discretization parameter N	Resolution of the regulatory interaction network: larger N captures finer-scale enhancer–promoter interactions and higher-order chromatin loops
Time T	Developmental or experimental time horizon; longer T corresponds to more accumulated regulatory history and greater epigenetic memory burden on the numerical scheme
Error decrease with N	Improving the resolution of the regulatory model consistently reduces the discrepancy between predicted and actual gene expression — the biological system is well-approximated by the mathematical framework
Error growth with T	Long developmental time horizons accumulate regulatory memory that is increasingly difficult to represent numerically — analogous to the growing complexity of the epigenome over developmental time

11.5. Time Effect — The Memory Burden of Developmental History

The most striking quantitative finding of the numerical analysis is the magnitude of temporal error amplification: from $T = 0.15$ to $T =0.85$, the maximum error grows by factors of 1,200 to 4,048 depending on the condition and kernel parameter. This is not a numerical artifact, it is a direct computational manifestation of a fundamental biological reality.

Through development, the epigenome accrues a record of increasing regulatory events. Each transcription factor binding event, each histone modification cycle, each DNA methylation change and each three-dimensional chromatin reorganisation contributes to the memory burden encoded in the Caputo fractional convolution kernel. The integral ${\int }_0^t{\left(t-\tau \right)}^{\alpha -1}f\left(\tau \right)d\tau$ remembers: it weights every past state by a power of the elapsed time, ensuring that the current transcriptional output of a gene reflects the accumulated influence of its entire regulatory history.

This is why differentiating a pluripotent stem cell into a terminally differentiated neuron is an irreversible process under normal conditions: the epigenetic memory accumulated over the course of differentiation, the growing t in the fractional integral, creates a regulatory state that is increasingly distant from the initial pluripotent condition and increasingly resistant to reprogramming. The numerical result that convergence rates fall from super linear values of 4–5 at $T=0.15$ to sub-linear values below 1.0 at $T=0.85$ is the mathematical expression of this progressive entrenchment of cell identity through epigenetic memory accumulation.

Critical implication for long-time gene regulation models: For gene regulatory simulations spanning long developmental timescales (corresponding to$T=0.45-0.85$in the dimensionless framework of Example 1), the numerical results demonstrate that temporal mesh refinement is at least as important as spatial mesh refinement. The marginal benefit of increasing $N$ from 6 to 12 diminishes as $T$ grows, while the error contributed by the fractional memory integration continues to accumulate. This implies that high-fidelity models of epigenetic inheritance across multiple cell divisions or developmental stages require adaptive temporal stepping strategies that allocate computational resources in proportion to the growing memory kernel integral

11.6. The Systematic Superiority of Nonlocal Conditions — A Biological Justification

The finding that nonlocal initial conditions produce smaller errors than local conditions in tested configurations, with reductions of 18–33%, carries a biological interpretation that goes beyond numerical convenience. It constitutes quantitative evidence for a fundamental principle of gene regulation: the current state of a gene cannot be fully explained by where it sits in the genome, but requires knowledge of how its regulatory state was established and maintained over developmental time.

Why nonlocal is biologically more accurate

Under local conditions, $w\left(x\right) = x⁴$ assigns the initial regulatory state of every gene as a function of chromosomal position alone. This corresponds to asking: given only the DNA sequence and chromosomal context of a locus at $t = 0$, what is its expression level? This is a well-posed but incomplete question. It captures position effects accurately but ignores the epigenetic marks, the$w(x,t)$ component, that have been laid down by the prior regulatory history of the cell. A gene at position x in a liver cell and the same gene at position x in a neuron are assigned identical initial conditions under the local model, yet they may have radically different expression levels because their chromatin states, their epigenetic memories, are entirely different.

The nonlocal condition $w(x,t)$ corrects this deficiency by encoding the spatiotemporal regulatory history of the locus directly into the initial condition. The parameter t in w(x,t) is not the future evolution variable; it is the compressed representation of the past: the record of all chromatin modifications, transcription factor binding events, and three-dimensional genome reorganizations that have shaped the current regulatory state of position x. When this history is incorporated into the initial condition, the numerical solution of the fractional integro-differential equation is more accurate, more convergent, and more physically faithful — precisely because the initial data is now compatible with the long-range, memory-laden structure of the fractional integral operator.

The connection to epigenetic inheritance

The nonlocal advantage grows with time. At $T=0.15$, the local condition is marginally better for ν = 0.02, but by $T=0.45$and $T=0.85$ the nonlocal condition consistently outperforms it for all kernel parameters. This temporal pattern mirrors the biology of epigenetic inheritance: at short developmental timescales, the positional component of gene regulation (local w(x)) dominates, because the epigenetic marks have not yet had time to diverge significantly between cells with different histories. At longer timescales, the accumulated regulatory history (nonlocal $w(x,t)$) becomes the primary determinant of gene expression, because the epigenome has been progressively shaped by the cell’s developmental experience. The crossover from local-competitive to nonlocal-superior behavior with increasing $T$ is therefore not merely a numerical observation, it is a computational demonstration of the developmental timescale at which epigenetic memory overtakes positional determinism as the governing force in gene regulation.

11.7. The Two Kernels as Two Modes of Genetic Regulation

The two kernel types are not alternatives to each other; they are complementary. In any real gene regulatory network, both cis and trans mechanisms operate simultaneously. The enhancer–promoter contacts of the Carleman type establish the spatial specificity of gene regulation which genes are regulated by which enhancers while the diffusion-mediated trans-regulatory interactions of the logarithmic type set the global transcriptional tone — the overall level of transcriptional activity in the nucleus. The fractional integro-differential framework of Example 1 accommodates both mechanisms within a single mathematical structure, and the numerical results confirm that the proposed method handles both kernel types with high accuracy and systematic convergence.

Table 14. Two Kernels as Two Modes of Genetic Regulation.

Property	Carleman kernel |y−τ|^(−β)	Logarithmic kernel ln|y−τ|
Singularity type	Algebraic (strong)	Logarithmic (weak)
Decay with distance	Power-law, rapid	Logarithmic, slow
Biological mechanism	Cis-regulatory: enhancer–promoter looping, TAD-mediated contact	Trans-regulatory: diffusion-mediated TF binding, nuclear condensate formation
Genomic scale	Kilobases to megabases	Genome-wide
Specificity	High (sequence-specific, position-dependent)	Low (motif-distributed, position-independent)
Convergence rate at short T	Very high (4–5 order)	High (5 order, Example 2)
Temporal sensitivity	High (rapid error growth with T)	Moderate (stable long-time performance)
Error at short T (N=12)	10⁻⁷ to 10⁻⁴	10⁻¹¹ (near machine precision)
Nonlocal advantage	Moderate to strong (18–33%)	Strong at medium T
Biological example	Hox gene clusters, β-globin locus control region, super-enhancers	Nuclear hormone receptors, pioneer transcription factors, phase-separated condensates

Table 14. Two Kernels as Two Modes of Genetic Regulation.

Property	Carleman kernel |y−τ|^(−β)	Logarithmic kernel ln|y−τ|
Singularity type	Algebraic (strong)	Logarithmic (weak)
Decay with distance	Power-law, rapid	Logarithmic, slow
Biological mechanism	Cis-regulatory: enhancer–promoter looping, TAD-mediated contact	Trans-regulatory: diffusion-mediated TF binding, nuclear condensate formation
Genomic scale	Kilobases to megabases	Genome-wide
Specificity	High (sequence-specific, position-dependent)	Low (motif-distributed, position-independent)
Convergence rate at short T	Very high (4–5 order)	High (5 order, Example 2)
Temporal sensitivity	High (rapid error growth with T)	Moderate (stable long-time performance)
Error at short T (N=12)	10⁻⁷ to 10⁻⁴	10⁻¹¹ (near machine precision)
Nonlocal advantage	Moderate to strong (18–33%)	Strong at medium T
Biological example	Hox gene clusters, β-globin locus control region, super-enhancers	Nuclear hormone receptors, pioneer transcription factors, phase-separated condensates

12. General Conclusion

The numerical results presented across nine convergence tables and six graphical figures establish, with comprehensive and systematic evidence, the accuracy, stability, and convergence of the proposed separation technique method combined with the Toeplitz Matrix Technique for the numerical solution of fractional integro-differential equations under both local and nonlocal initial conditions.

Three principal findings emerge from the data. First, the method is globally convergent under all tested conditions: increasing the discretization parameter N from 3 to 6 to 12 reduces the maximum absolute error without exception across all nine tables, for both values of ν, at all three time levels, and under both local and nonlocal conditions. This universality of convergence is the most fundamental validation that the TMT is a reliable and robust numerical scheme for this class of equations.

Second, time is the dominant determinant of numerical accuracy. In Example 1. The transition from$T=0.15$ to $T=0.85$ increases the maximum error by two to six orders of magnitude depending on the example and configuration, while the convergence rate declines from superlinear values of 4–5 at short times to sub-linear values below 1.0 at long times. This behavior is a direct consequence of the cumulative memory effect of the Caputo fractional derivative, whose power-law integration kernel accumulates all past regulatory history. For biological applications involving long developmental timescales or persistent epigenetic states, this finding implies that high-accuracy fractional simulations require temporal mesh refinement strategies that adapt to the growing memory burden, not merely spatial refinement.

Third, nonlocal initial conditions yield systematically superior numerical accuracy compared to local conditions, with error reductions of 20–35% observed across comparable configurations. This finding reflects a deep compatibility between the nonlocal structure of the initial data and the global basis representation of the Toeplitz method: when the initial condition encodes the same long-range spatiotemporal dependencies as the integral operator, the projection residual is smaller and the convergence is faster. From a biological standpoint, this result validates the use of epigenetically informed initial conditions — conditions that carry the regulatory history of the cell as the appropriate starting point for fractional models of gene expression dynamics.

The comparison of Carleman and logarithmic kernels reveals that both kernel types are faithfully and accurately handled by the proposed scheme, but with complementary performance profiles. The Carleman kernel, with its stronger algebraic singularity, produces higher convergence rates at short times (reaching fifth order in Example 2, Table 7) but greater temporal sensitivity. The logarithmic kernel, with its weaker singularity, produces more moderate but stable convergence rates across the full temporal range. In the context of gene regulation, these two kernels correspond respectively to cis-regulatory enhancer-promoter contacts (strong, distance-dependent, algebraically decaying) and trans-regulatory diffusion-mediated interactions (weak, spatially diffuse, logarithmically decaying). The ability of the proposed method to handle both kernel types with high accuracy positions it as a versatile computational tool for modeling the full spectrum of gene regulatory mechanisms.

The surface plots in Figure 1 through 4 demonstrate that the approximate solutions are smooth, structurally stable, and physically interpretable under a range of nonlocal correction functions. The qualitative contrast between the monotone decay of Figure 1, the sign-changing saddle structure of Figure 2, the antisymmetric growth of Figure 3, and the monotone ramp of Figure 4 reveals that the nonlocal correction function is the primary determinant of the global solution morphology, a finding that parallels the biological observation that the epigenetic state of a gene's regulatory landscape, rather than its coding sequence, determines the spatial and temporal pattern of its expression. In Figure 5, the form of potential function in local condition for logarithmic kernel. The two-dimensional error profiles of confirm that the spatial error is concentrated near the boundary of the approximation interval and remains within the theoretically predicted bounds at $N = 8$. The near-perfect visual coincidence of Approx1 and Approx2 across the full spatial domain at fixed $T = 0.45$ provides a strong graphical confirmation of solution reproducibility and numerical stability.

In summary, the proposed fractional computational framework, validated herein across a comprehensive parameter space, constitutes a mathematically rigorous and physically justified tool for modeling genetic gene regulation. It captures simultaneously the fast transcriptional dynamics and slow epigenetic memory through its dual-order fractional structure, the spatial specificity of cis-regulatory interactions through the Carleman kernel, the diffuse long-range action of trans-regulatory factors through the logarithmic kernel, the positional determinism of classical genetics through local initial conditions, and the history-dependence of epigenetic inheritance through nonlocal initial conditions. The systematic superiority of the nonlocal framework at all tested configurations constitutes quantitative evidence that epigenetic memory is not merely a biological hypothesis but a computationally measurable and numerically significant component of gene regulatory dynamics.