Anomalous diffusion is a well-documented phenomenon in biological tissues, where complex microstructural features such as cellular crowding, membrane barriers, and heterogeneous extracellular matrices give rise to non-Gaussian transport and long-range memory effects. Classical diffusion models and constant-order fractional equations are often insufficient to capture the spatial variability of diffusion mechanisms observed across different tissue regions. Variable-order fractional diffusion models provide a natural and physiologically meaningful framework to describe such heterogeneous transport processes.

In this work, we propose an efficient and energy-preserving numerical method for a class of time–space fractional diffusion equations with spatially variable fractional order, motivated by anomalous transport in heterogeneous biological tissues. The model incorporates a Caputo time-fractional derivative and a piecewise-defined spatial fractional diffusion operator, enabling the coexistence of distinct diffusion regimes associated with local microstructural properties.

The numerical scheme combines a memory-efficient temporal discretisation of the fractional derivative with a stable spatial approximation of the variable-order operator. A discrete energy functional is constructed, and unconditional stability of the fully discrete scheme is rigorously established through preservation of the dissipative energy structure of the continuous problem. Convergence is proven under mild regularity assumptions.

Numerical experiments in two-dimensional heterogeneous tissue-like domains validate the theoretical results and demonstrate biologically relevant diffusion phenomena, including spatially dependent propagation speeds and effective diffusion barriers induced by microstructural heterogeneity. These results illustrate the limitations of constant-order models and highlight the relevance of variable-order fractional formulations for quantitative modelling of anomalous transport in complex biological tissues.