Efficient and Structure-Preserving Numerical Methods for Time--Space Fractional Diffusion in Heterogeneous Biological Tissues

1. Introduction

1.1. Biomedical Motivation

Diffusion processes play a central role in numerous physiological and pathological mechanisms, including nutrient transport, metabolite exchange, signalling molecule propagation, and drug delivery within biological tissues. Classical diffusion models, derived from Fick’s laws, rely on assumptions of Brownian motion, spatial homogeneity, and the absence of long-range temporal correlations. Although these assumptions are mathematically convenient, they often fail to capture transport phenomena observed in living tissues [1,2].

Biological tissues are characterised by a highly complex and heterogeneous microstructure, including cellular membranes, organelles, extracellular matrix fibres, vascular networks, and varying pore geometries. These structural features introduce obstacles, trapping effects, transient binding, and spatial heterogeneity, all of which strongly influence particle motion. As a consequence, experimental measurements frequently reveal deviations from the linear mean squared displacement predicted by classical Fickian diffusion [3,4].

Anomalous diffusion has been reported in a wide range of biological systems using experimental techniques such as fluorescence recovery after photobleaching (FRAP), diffusion-weighted magnetic resonance imaging (DW–MRI), and single-particle tracking. Subdiffusive behaviour is commonly observed in crowded intracellular environments and soft biological matter, where viscoelasticity and molecular crowding dominate transport dynamics [5,6]. In contrast, superdiffusive transport may emerge in tissues exhibiting active processes, directed motion, or long-range structural correlations [7].

Fractional-order diffusion models provide a mathematically consistent and physically interpretable framework for describing such non-Fickian transport phenomena. Temporal fractional derivatives naturally incorporate memory effects and long-time correlations, while spatial fractional operators account for nonlocal transport and long-range interactions induced by tissue heterogeneity [8,9]. Importantly, the fractional orders appearing in these models can often be linked to experimentally measurable properties, establishing a meaningful connection between mathematical formulation and biological interpretation [10].

From an engineering and clinical perspective, accurate modelling of anomalous diffusion is essential for improving predictive simulations at the tissue scale, optimising drug delivery strategies, and supporting the design of biomedical devices and therapeutic protocols. These considerations provide a strong motivation for the development of reliable, efficient, and mathematically rigorous computational methods for fractional diffusion equations arising in biological tissue modelling.

1.2. Mathematical Background

Fractional calculus extends classical differential operators by allowing derivatives and integrals of non-integer order. Unlike standard integer-order derivatives, fractional derivatives are inherently nonlocal operators, as they involve integration over an interval in time or space. This nonlocality provides a natural mathematical mechanism to incorporate memory effects and long-range interactions, which are essential features of anomalous transport processes observed in complex media [8,9].

In the temporal domain, fractional derivatives account for history-dependent dynamics. For instance, the Caputo fractional derivative of order $0<\alpha <1$ introduces a convolution kernel with algebraic decay, reflecting long-term memory effects and power-law waiting times between particle movements. Such behaviour is a hallmark of subdiffusive transport and cannot be captured by classical first-order time derivatives [3,11]. From a modelling perspective, time-fractional diffusion equations arise naturally as macroscopic limits of continuous-time random walk processes with heavy-tailed waiting-time distributions.

In the spatial domain, fractional operators such as the fractional Laplacian or Riesz derivative generalise the classical second-order Laplace operator by incorporating nonlocal spatial interactions. These operators describe transport processes characterised by long jumps or Lévy flights, which are particularly relevant in heterogeneous or porous media. In biological tissues, spatial fractional derivatives provide a mathematical description of transport influenced by structural irregularities, connectivity across scales, and long-range correlations induced by the underlying microarchitecture [12,13].

Throughout this work, we adopt the spectral definition of the fractional Laplacian ${(-\Delta )}^{\beta /2}$, $\beta \in (1,2]$, defined via the eigenpairs of the classical Laplacian subject to the prescribed boundary conditions. This choice ensures a rigorous variational formulation in fractional Sobolev spaces and is well suited for finite element discretisations on bounded domains.

Time–space fractional diffusion equations combine these two generalisations into a unified mathematical framework, leading to models of the form

$${{}_0{}^{C}D}_t^{\alpha }u(x,t)={\mathcal{L}}^{\beta }u(x,t)+f(x,t),$$

where ${{}_0{}^{C}D}_t^{\alpha }$ denotes a fractional temporal derivative and ${\mathcal{L}}^{\beta }$ represents a fractional spatial operator of order $\beta$. Classical diffusion equations are recovered as special cases when $\alpha =1$ and $\beta =2$, highlighting the role of fractional models as genuine extensions rather than ad hoc modifications of standard theory [14].

From a mathematical standpoint, fractional diffusion equations pose significant analytical and computational challenges. The presence of nonlocal operators complicates the analysis of well-posedness, regularity, and stability, while also increasing the computational cost of numerical approximations. Nevertheless, these models offer remarkable flexibility: the fractional orders $\alpha$ and $\beta$ act as tunable parameters that encode medium heterogeneity, memory strength, and transport modality, enabling precise adaptation to experimentally observed phenomena [10,15].

In the context of biological tissue modelling, the mathematical structure of time–space fractional diffusion equations provides a principled way to bridge microscopic transport mechanisms and macroscopic tissue-scale behaviour. This dual interpretability, combined with their strong theoretical foundations, makes fractional diffusion models particularly attractive for the quantitative description of anomalous transport in heterogeneous biological media and motivates the development of robust analytical and numerical tools for their solution.

1.3. Related Work

Over the past two decades, fractional diffusion equations have attracted considerable attention, leading to the development of a wide variety of analytical and numerical methods. Early numerical approaches largely focused on finite difference schemes for time-fractional diffusion equations, where standard discretisations of the Caputo derivative were combined with classical spatial operators. While these methods established a foundation for numerical approximation, they typically exhibited high computational and memory costs due to the global-in-time nature of fractional derivatives [8,11].

To address these limitations, several authors proposed improved finite difference and finite element methods for time–space fractional diffusion equations. Implicit and semi-implicit schemes have been developed to enhance numerical stability, often accompanied by convergence analyses under restrictive regularity assumptions [16,17]. Despite these advances, many schemes remain computationally expensive, as the memory requirement grows linearly with the number of time steps, making large-scale or long-time simulations impractical.

Spatial discretisation of fractional operators has been investigated using a variety of techniques, including Grünwald–Letnikov approximations, matrix transfer techniques, finite element methods, and spectral approaches. Finite element formulations offer flexibility in complex geometries but typically involve dense stiffness matrices when fractional Laplacians are employed, leading to significant computational overhead [18,19]. Spectral methods provide high accuracy for smooth solutions, yet their applicability is limited in heterogeneous biological domains with irregular boundaries [20].

More recently, structure-preserving numerical methods have gained increasing attention. Such approaches aim to maintain qualitative properties of the continuous model, including positivity, monotonicity, and mass conservation. However, many existing fractional schemes focus primarily on accuracy and convergence, without ensuring the preservation of physically meaningful invariants, which are particularly important in biomedical applications [21].

To mitigate the high memory cost inherent to fractional derivatives, fast algorithms based on sum-of-exponentials approximations, short-memory principles, and adaptive time stepping have been introduced. These methods significantly reduce storage and computational complexity but often complicate stability analysis or introduce additional approximation errors that are not fully quantified [22,23].

Despite the substantial progress in numerical methods for fractional diffusion equations, several challenges remain open. There is a clear need for numerical schemes that simultaneously achieve memory efficiency, unconditional stability, rigorous convergence guarantees, and preservation of key physical structures. Furthermore, relatively few studies explicitly address the demands imposed by heterogeneous biological tissues, where variability in material properties and complex boundary conditions exacerbate numerical difficulties. These gaps in the existing literature motivate the development of the computational framework proposed in this work.

1.4. Contributions and Novelty

The main contributions of this work can be summarised as follows:

• the development of a fully discrete implicit numerical scheme for time–space fractional diffusion equations with heterogeneous coefficients;
• the proof of unconditional stability with respect to the time step, fractional order $\alpha$, and spatial discretisation parameter h;
• a reduction of the memory complexity of the classical L1 scheme from $\mathcal{O}\left(N\right)$ to $\mathcal{O}(logN)$ via a sum-of-exponentials approximation;
• numerical validation of the method in heterogeneous media with discontinuous diffusion coefficients, relevant to biological tissue modelling.

2. Mathematical Model of Anomalous Diffusion in Biological Tissues

2.1. Governing Equation

Let $\Omega \subset {\mathbb{R}}^d$, with $d\ge 1$, be a bounded open domain representing a biological tissue region, and let $T>0$ denote the final observation time. We consider the following time–space fractional diffusion equation:

$${{}_0{}^{C}D}_t^{\alpha }u(x,t)+{(-\nabla ·(\kappa \left(x\right)\nabla ))}^{\beta /2}u(x,t)=f(x,t),(x,t)\in \Omega \times (0,T].$$

(1)

where $u(x,t)$ denotes the concentration of a diffusing quantity, such as a solute or tracer, within the tissue. Here, ${(-\nabla ·(\kappa \left(x\right)\nabla ))}^{\beta /2}$ denotes a variable-coefficient fractional diffusion operator defined in the spectral sense, which reduces to the classical diffusion operator when $\beta =2$. The operator ${{}_0{}^{C}D}_t^{\alpha }$ is the Caputo fractional derivative of order $0<\alpha <1$, defined by

$${{}_0{}^{C}D}_t^{\alpha }u\left(t\right)={\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\int }_0^t{(t-s)}^{-\alpha }{\displaystyle \frac{\partial u\left(s\right)}{\partial s}}ds,$$

(2)

which naturally incorporates memory effects via a convolution kernel with algebraic decay [8,9].

The spatial operator ${\nabla }^{\beta }$ denotes a space-fractional diffusion operator of order $1<\beta \le 2$, typically interpreted as the Riesz fractional gradient or, equivalently, as a fractional Laplacian in divergence form. This operator accounts for nonlocal spatial interactions and is particularly suitable for modelling transport in heterogeneous and structurally complex media [12,19]. The coefficient $\kappa \left(x\right)$ represents a spatially varying diffusivity tensor that reflects heterogeneity in tissue properties. The source term $f(x,t)$ models production, consumption, or external input of the diffusing species [3,10].

Equation (??) generalises the classical diffusion equation. In the limiting case $\alpha =1$ and $\beta =2$, the standard Fickian diffusion model is recovered, highlighting the fractional formulation as a mathematically consistent extension rather than a phenomenological modification [14].

2.2. Biological Interpretation

The fractional temporal order $\alpha$ reflects the strength of memory effects in the transport process. Values of $\alpha <1$ correspond to subdiffusive behaviour, which is commonly observed in biological tissues due to molecular crowding, binding interactions, and viscoelasticity of the intracellular and extracellular environments [5,6]. In this context, the Caputo derivative provides a causal representation of memory, making it suitable for modelling physiological processes with well-defined initial states.

The spatial fractional order $\beta$ characterises the degree of nonlocal spatial transport. Values $\beta <2$ indicate the presence of long-range particle displacements or Lévy-type transport mechanisms, which may arise from heterogeneous tissue architecture, porous structures, or multiscale connectivity within the biological medium [4,7]. The combination of time- and space-fractional operators allows the model to capture a wide spectrum of anomalous diffusion behaviours observed experimentally.

The spatially varying diffusivity $\kappa \left(x\right)$ further enhances the biological realism of the model by accounting for tissue heterogeneity, such as variations in cellular density, extracellular matrix composition, or perfusion. Together, the parameters $(\alpha ,\beta ,\kappa )$ form a compact yet expressive representation of tissue-scale transport dynamics [10,15].

2.3. Initial and Boundary Conditions

To complete the mathematical model, equation (??) is supplemented with an initial condition

$$u(x,0)=u_0\left(x\right),x\in \Omega ,$$

(3)

where $u_0\left(x\right)$ represents the initial distribution of the diffusing quantity, assumed to be non-negative and sufficiently regular.

Boundary conditions are chosen to reflect physiological constraints and depend on the specific biomedical application. Common choices include homogeneous Dirichlet conditions,

$$u(x,t)=0,x\in \partial \Omega ,t>0,$$

(4)

which may model complete absorption or clearance at tissue boundaries, and homogeneous Neumann-type conditions,

$$\left(\kappa \left(x\right){\nabla }^{\beta }u(x,t)\right)·n\left(x\right)=0,x\in \partial \Omega ,t>0,$$

(5)

representing no-flux or impermeable boundaries [13]. We remark that, for nonlocal fractional operators, classical Neumann boundary conditions require careful interpretation. In the present work, boundary conditions are understood in the variational sense associated with the chosen spectral definition of the operator.

More general boundary conditions, including mixed or nonlocal formulations, can also be considered to describe exchange with surrounding tissues or vascular compartments. The well-posedness and numerical treatment of such conditions in the fractional setting require particular care due to the nonlocal nature of the spatial operator [18,19]. These modelling choices are discussed in greater detail in the numerical experiments section.

3. Functional Setting and Analytical Properties

3.1. Fractional Sobolev Spaces

The analysis of time–space fractional diffusion equations requires an appropriate functional framework that accounts for the nonlocal nature of the involved operators. Let $\Omega \subset {\mathbb{R}}^d$ be a bounded Lipschitz domain. For $s\in (0,1)$ and $1\le p<\infty$, the fractional Sobolev space $W^{s,p}(\Omega )$ is defined as

$$W^{s,p}(\Omega )=\left\{v\in L^p(\Omega ):{\int }_{\Omega }{\int }_{\Omega }{\displaystyle \frac{{|v\left(x\right)-v\left(y\right)|}^p}{{|x-y|}^{d+sp}}}dxdy<\infty \right\},$$

(6)

endowed with the norm

$${\parallel v\parallel }_{W^{s,p}(\Omega )}={\left({\parallel v\parallel }_{L^p(\Omega )}^p+{\int }_{\Omega }{\int }_{\Omega }{\displaystyle \frac{{|v\left(x\right)-v\left(y\right)|}^p}{{|x-y|}^{d+sp}}}dxdy\right)}^{1/p}.$$

(7)

In the Hilbertian case $p=2$, we write $H^s(\Omega ):=W^{s,2}(\Omega )$.

Fractional Sobolev spaces provide the natural variational setting for space-fractional operators, such as the fractional Laplacian. In particular, for $\beta \in (1,2]$, the operator ${\nabla }^{\beta }$ is typically associated with the space $H^{\beta /2}(\Omega )$, where weak formulations and energy estimates can be rigorously established [19,24]. These spaces generalise classical Sobolev spaces and retain key properties such as completeness, reflexivity, and compact embeddings under suitable conditions.

In the temporal dimension, the fractional Caputo derivative naturally leads to solution spaces involving $L^2(0,T;H^{\beta /2}(\Omega ))$ together with weighted Sobolev spaces in time that reflect the history dependence of the solution. Detailed discussions of fractional Sobolev regularity in time can be found in [11,25].

3.2. Well-Posedness

We now address the well-posedness of the time–space fractional diffusion problem introduced in Section 2. Let $u_0\in L^2(\Omega )$, $f\in L^2(0,T;L^2(\Omega ))$, and assume that the diffusion coefficient $\kappa \left(x\right)$ satisfies

$$0<{\kappa }_{min}\le \kappa \left(x\right)\le {\kappa }_{max}<\infty \mathrm{for}\mathrm{almost}\mathrm{every}x\in \Omega .$$

Definition 1. 

A function $u\in L^2(0,T;H^{\beta /2}(\Omega ))$ is called a weak solution of the time–space fractional diffusion problem if $u(·,0)=u_0$ and, for all test functions $v\in H^{\beta /2}(\Omega )$ and almost every $t\in (0,T)$,

$$\left({{}_0{}^{C}D}_t^{\alpha }u\left(t\right),v\right)+a(u\left(t\right),v)=(f\left(t\right),v),$$

(8)

where $(·,·)$ denotes the $L^2(\Omega )$ inner product and

$$a(u,v):={\int }_{\Omega }\kappa \left(x\right){(-\Delta )}^{\beta /4}u{(-\Delta )}^{\beta /4}vdx.$$

This bilinear form is symmetric, coercive and continuous on $H^{\beta /2}(\Omega )$ under the assumption $0<{\kappa }_{min}\le \kappa \left(x\right)\le {\kappa }_{max}<\infty$.

Theorem 1. 

Under the above assumptions on f, $u_0$, and $\kappa \left(x\right)$, the time–space fractional diffusion problem admits a unique weak solution

$$u\in L^2(0,T;H^{\beta /2}(\Omega )),$$

which depends continuously on the data.

Proof 

(Sketch of proof). The proof follows from standard arguments combining fractional energy estimates with the theory of maximal monotone operators. The coercivity and boundedness of the bilinear form $a(·,·)$ in $H^{\beta /2}(\Omega )$ ensure spatial well-posedness. The properties of the Caputo derivative yield an appropriate fractional Grönwall inequality, which guarantees uniqueness and continuous dependence on the data [23,25]. Existence is obtained via Galerkin approximations and compactness arguments adapted to the fractional setting [14]. □

This well-posedness result establishes a rigorous analytical foundation for the numerical methods developed in the subsequent sections.

Throughout the numerical analysis, we assume that the exact solution u of problem (1) satisfies

$$u\in C([0,T];H^{\beta /2}(\Omega )),{\partial }_{t}u\in L^1(0,T;H^{\beta /2}(\Omega )),$$

and there exists a constant $C>0$ such that

$$\parallel {\partial }_t{u\left(t\right)\parallel }_{H^{\beta /2}(\Omega )}\le C{t}^{\alpha -1},t\in (0,T].$$

This weak singularity at the initial time is a characteristic feature of time-fractional diffusion equations and is essential for deriving sharp temporal error estimates.

4. Numerical Method

This section presents a fully discrete numerical method for the time–space fractional diffusion problem (??). Particular attention is paid to the nonlocal nature of both the temporal and spatial operators, with the goal of constructing a scheme that is stable, convergent, and memory efficient, while remaining suitable for heterogeneous biological tissues.

4.1. Temporal Discretisation

Let $0=t_0<t_1<\dots <t_N=T$ be a uniform partition of the time interval $[0,T]$ with time step $\Delta t=T/N$, and denote $u^n\left(x\right)\approx u(x,t_n)$. The Caputo fractional derivative at time $t=t_n$ is defined by

$${{}_0{}^{C}D}_t^{\alpha }u(x,t_n)={\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\int }_0^{t_n}{(t_n-s)}^{-\alpha }{\displaystyle \frac{\partial u(x,s)}{\partial s}}ds.$$

A classical and robust approximation of this operator is provided by the L1 scheme, which is derived by approximating ${\partial }_{t}u$ by a backward finite difference on each subinterval $[t_{k-1},t_k]$. This yields

$${{}_0{}^{C}D}_t^{\alpha }u(x,t_n)\approx {\displaystyle \frac{1}{\Gamma (2-\alpha ){(\Delta t)}^{\alpha }}}\sum _{k=0}^{n-1}{b}_k\left(u^{n-k}\left(x\right)-u^{n-k-1}\left(x\right)\right),$$

(9)

where the convolution weights are given by

$$b_k={(k+1)}^{1-\alpha }-k^{1-\alpha },k\ge 0.$$

The $L^1$ scheme is first-order accurate in time for sufficiently smooth solutions and has the advantage of unconditional stability when combined with implicit spatial discretisations. However, a well-known drawback is its $\mathcal{O}\left(N\right)$ memory requirement per spatial degree of freedom, which becomes prohibitive for long-time simulations typical in biological diffusion problems.

To reduce memory and computational costs, we adopt a memory-efficient convolution strategy based on truncation and compression of the history term. Specifically, the convolution sum in (9) is decomposed into a recent-history part, computed exactly, and a distant-history part, approximated using a sum-of-exponentials (SOE) representation of the kernel $t^{-\alpha }$. This approach preserves the accuracy and stability of the $L^1$ scheme while reducing the memory complexity from $\mathcal{O}\left(N\right)$ to $\mathcal{O}(logN)$, which is essential for realistic biomedical simulations over long time horizons.

4.2. Spatial Discretisation

We now describe the discretisation of the space-fractional diffusion operator. For clarity of exposition, we focus on the one-dimensional case $\Omega =(a,b)$; higher-dimensional extensions follow analogously.

The fractional diffusion operator is interpreted as a Riesz-type fractional Laplacian of order $\beta \in (1,2]$, written in weak form as

$${\mathcal{L}}^{\beta }u=-\nabla ·\left(\kappa \left(x\right){\nabla }^{\beta -1}u\right).$$

We discretise this operator using a finite element framework, which is particularly well suited to heterogeneous tissues and complex geometries. Let ${\mathcal{T}}_h$ be a conforming triangulation of $\Omega$ with mesh size h, and let $V_h\subset H^{\beta /2}(\Omega )$ be a space of continuous, piecewise linear finite elements satisfying the prescribed boundary conditions.

The semi-discrete weak formulation reads: find $u_h\left(t\right)\in V_h$ such that, for all $v_h\in V_h$,

$$\left({{}_0{}^{C}D}_t^{\alpha }{u}_h\left(t\right),v_h\right)+a(u_h\left(t\right),v_h)=(f\left(t\right),v_h),$$

(10)

where the bilinear form

$$a(u,v)={\int }_{\Omega }\kappa \left(x\right){\nabla }^{\beta /2}u·{\nabla }^{\beta /2}vdx$$

is coercive and bounded on $H^{\beta /2}(\Omega )$.

The resulting stiffness matrix is dense due to the nonlocal nature of the fractional operator. However, its structure is symmetric positive definite, and efficient matrix compression techniques or truncated interaction strategies can be employed to mitigate computational cost without destroying stability.

4.3. Fully Discrete Scheme

Combining the temporal discretisation (9) with the finite element spatial approximation leads to the following fully discrete scheme.

For $n=1,2,\dots ,N$, find $u_h^n\in V_h$ such that, for all $v_h\in V_h$,

$${\displaystyle \frac{1}{\Gamma (2-\alpha ){(\Delta t)}^{\alpha }}}\sum _{k=0}^{n-1}{b}_k\left(u_h^{n-k}-u_h^{n-k-1},v_h\right)+a(u_h^n,v_h)=(f^n,v_h).$$

(11)

This scheme is implicit with respect to the spatial operator, which is essential to ensure unconditional stability, particularly for small values of $\alpha$. The resulting linear system at each time step involves the same stiffness matrix, allowing for efficient reuse of factorisations or preconditioners.

4.4. Computational Complexity

We conclude this section by analysing the computational cost of the proposed method.

Let M denote the number of spatial degrees of freedom and N the number of time steps. A direct implementation of the L1 scheme would require $\mathcal{O}\left(MN\right)$ storage and $\mathcal{O}\left(MN\right)$ operations. By employing a compressed SOE-based history approximation, the memory requirement is reduced to $\mathcal{O}(MlogN)$ and the total computational cost to $\mathcal{O}(MNlogN)$.

Importantly, the spatial discretisation cost is dominated by the solution of a linear system with a dense but structured matrix. In practice, this cost can be significantly reduced using hierarchical matrices or fast multipole–type techniques, which are particularly effective for fractional operators.

Overall, the proposed numerical method achieves a favourable balance between accuracy, stability, and computational efficiency, making it suitable for large-scale simulations of anomalous diffusion in heterogeneous biological tissues.

4.5. Stability of the Fully Discrete Scheme

Theorem 2 

(Unconditional stability in $L^2(\Omega )$). Let $u_h^n\in V_h$ be the solution of the fully discrete scheme (11). Under the assumptions on the diffusion coefficient $\kappa \left(x\right)$, the scheme is unconditionally stable in the sense that

$$\parallel u_h^n{\parallel }_{L^2(\Omega )}^2+\Delta t^{\alpha }\sum _{k=1}^n{\parallel u_h^k\parallel }_{H^{\beta /2}(\Omega )}^2\le C\left(\parallel u_0{\parallel }_{L^2(\Omega )}^2+\underset{1\le k\le n}{max}{\parallel f^k\parallel }_{L^2(\Omega )}^2\right),$$

where the constant $C>0$ is independent of h, $\Delta t$, and n.

Proof 

(Sketch of proof). The result follows from the positivity and monotonicity of the L1 convolution weights, the coercivity of the bilinear form $a(·,·)$, and a discrete fractional Grönwall inequality. □

4.6. Error Analysis

4.6.1. Temporal Discretisation Error

Theorem 3 

(Temporal convergence of the L1 scheme). Assume that the spatial discretisation is exact and that the regularity assumption holds. Then the L1 approximation of the Caputo derivative satisfies

$$\parallel u\left(t_n\right)-u^n{\parallel }_{L^2(\Omega )}\le C\Delta t^{min(1,2-\alpha )},1\le n\le N,$$

where the constant $C>0$ is independent of the time step $\Delta t$.

This rate reflects the reduced temporal regularity induced by the fractional time derivative and is optimal under the assumed solution regularity.

4.6.2. Spatial Discretisation Error

Theorem 4 

(Spatial convergence). Let $u\left(t\right)\in H^{\beta /2+s}(\Omega )$ for some $0<s\le 1-\beta /2$. Then the finite element approximation satisfies

$$\parallel u\left(t\right)-u_h{\left(t\right)\parallel }_{L^2(\Omega )}\le C{h}^s{\parallel u\left(t\right)\parallel }_{H^{\beta /2+s}(\Omega )},$$

where the constant $C>0$ is independent of the mesh size h.

4.6.3. Fully Discrete Error with SOE Approximation

Theorem 5 

(Fully discrete error estimate). Assume that the kernel $t^{-\alpha }$ is approximated by a sum-of-exponentials representation with accuracy ${\epsilon }_{\mathrm{SOE}}$. Then the fully discrete solution $u_h^n$ satisfies

$$\parallel u\left(t_n\right)-u_h^n{\parallel }_{L^2(\Omega )}\le C\left(\Delta t^{min(1,2-\alpha )}+h^s+{\epsilon }_{\mathrm{SOE}}\right),1\le n\le N,$$

where $C>0$ is independent of h, $\Delta t$, and n.

This estimate explicitly quantifies the trade-off between temporal accuracy, spatial resolution, and memory compression.

5. Numerical Experiments

In all experiments, the numerical solution remained bounded for arbitrarily small time step sizes, providing numerical evidence of the unconditional stability established in Theorem 2.

5.1. Temporal Convergence

We first assess the temporal convergence properties of the proposed numerical method. It is worth noting that the convergence rates observed in this experiment depend on the temporal regularity of the exact solution. Although the theoretical analysis provides a lower bound on the convergence rate, higher orders may be observed in practice for sufficiently smooth solutions. We consider the one-dimensional domain $\Omega =(0,1)$ and prescribe the manufactured solution

$$u(x,t)=t^{\alpha }sin\left(\pi x\right),$$

for which the corresponding source term $f(x,t)$ is computed analytically. A sufficiently fine spatial mesh is employed so that spatial discretisation errors are negligible compared to temporal errors. Table 1 reports the $L^2$-errors at the final time T for decreasing time step sizes $\Delta t$.

The temporal convergence behaviour of the proposed L1 scheme is illustrated in Figure 1, which displays the $L^2$-error at the final time as a function of the time step size $\Delta t$ on a log–log scale.

5.2. Spatial Convergence

To investigate spatial convergence, we fix a sufficiently small time step and successively refine the spatial mesh. The experiment is repeated for different values of the fractional diffusion order $\beta$. This strategy ensures that the reported errors predominantly reflect the spatial discretisation.

The spatial convergence behaviour of the proposed finite element discretisation is shown in Figure 2 for several values of $\beta$.

The observed convergence rates remain below first order, which is consistent with the limited spatial regularity of fractional diffusion problems.

5.3. Efficiency and Memory Reduction

We next compare the classical $L^1$ scheme with the proposed sum-of-exponentials (SOE) based memory compression strategy. For increasing final times T and a fixed time step size $\Delta t$, we record the memory required to store the temporal history. This experiment is designed to assess the asymptotic memory complexity of the two approaches.

Figure 3 illustrates the memory usage as a function of the number of time steps. While the classical $L^1$ scheme requires storage of the entire solution history and therefore exhibits linear growth in memory, the SOE-based strategy achieves a logarithmic memory footprint.

Additional tests (not shown) confirmed that the use of the SOE approximation does not introduce any noticeable degradation in accuracy compared to the classical $L^1$ scheme, provided that the SOE tolerance is chosen sufficiently small.

5.4. Heterogeneous Diffusion

Finally, we investigate anomalous diffusion in a heterogeneous medium characterised by a discontinuous diffusion coefficient

$$\kappa \left(x\right)=\left\{\begin{array}{cc}1,\hfill & x<0.5,\hfill \\ 10,\hfill & x\ge 0.5.\hfill \end{array}\right.$$

This configuration models a sharp transition in material properties and is representative of heterogeneous biological tissues. The effect of spatial heterogeneity on the diffusion process is illustrated in Figure 4, which shows solution profiles at different time instances.

These results directly support the modelling motivation discussed in Section 1.1, highlighting the ability of fractional diffusion models to capture heterogeneous transport mechanisms.

6. Conclusion

In this work, we developed and analysed an efficient and structure-preserving numerical framework for time–space fractional diffusion equations arising in heterogeneous biological tissues. The proposed approach combines an implicit finite element discretisation of the fractional spatial operator with a memory-efficient L1 approximation of the Caputo time derivative, yielding a stable and robust scheme suitable for long-time simulations.

Rigorous stability and convergence analyses were established, including a fully discrete error estimate that explicitly quantifies the contributions of temporal discretisation, spatial approximation, and kernel compression. In particular, the theoretical results demonstrate that the incorporation of a sum-of-exponentials approximation substantially reduces memory requirements while preserving numerical accuracy, thereby achieving a favourable balance between efficiency and reliability.

Comprehensive numerical experiments confirmed the predicted convergence rates and provided numerical evidence supporting the theoretical findings. Moreover, simulations in heterogeneous media illustrated the strong interplay between spatial variability and fractional transport mechanisms, directly supporting the modelling motivation presented in Section 1.1 and highlighting the limitations of classical diffusion models in complex biological tissues.

Overall, the proposed framework constitutes a mathematically sound and computationally scalable tool for the simulation of anomalous diffusion in heterogeneous media. Future work will focus on adaptive time-stepping strategies, higher-order temporal discretisations, and extensions to three-dimensional biomedical geometries informed by experimental and imaging data.