## 1 Introduction  
Because of advancements in microwave, laser, focused ultrasound, and radiofrequency technologies, many modern thermo-therapeutics are now widely used in clinical treatment. For example, in thermal therapy, an objective lens focuses a laser on a tumor. One of the most difficult problems in thermal therapy is delivering the right amount of heat energy to the diseased tissue while avoiding damage to healthy tissue. As a result, there is a pressing need to comprehend how temperature/stress fields influence the kinetics of thermal therapy. [PERSON] and [PERSON] [1] demonstrated that deformation caused by heating and cooling can cause pain sensation. Accurate predictions of heat and mechanical reactions, as well as thermal damage in biological tissue, are thus essential for treatment planning and the development of novel therapeutic heating systems. Because of the inherent properties of biological tissue, heat transfer analysis in living biological tissue is a complex physiological procedure. Heat transfer analysis in living biological tissue is a difficult physiological process due to biological tissue's inherent properties, such as blood circulation, sweating, metabolic heat generation, and heat dissipation via hair or fur. [PERSON] [2] used his bioheat transfer model to simulate the temperature profile in the human forearm to adequately characterize this complex event. Other bioheat transfer models were developed to address the shortcomings of [PERSON] model, which does not consider the blood velocity field or the geometry of blood vessels [3; 4; 5; 6; 7]. Although these models are more detailed and precise than [PERSON]', their complexities make them difficult to apply in practice. On the other hand, because [PERSON] model is simple and has a small number of material parameters, it has piqued the interest of many researchers [8]. Because of its non-homogeneous structure and composition, soft tissue's mechanical response in vivo is heterogeneous, anisotropic, nonlinear, and viscoelastic. It is influenced by a variety of factors, including age, gender, location, and hydration. [PERSON] [9] and [PERSON] [10] proposed a single-phase lag (SPL) related to heat flow to address a limitation of the linear bioheat transfer model. The thermal wave bioheat transfer (TWBHT) model was obtained by combining the SPL model with the energy equation. [PERSON] [11,12] extended the SPL model by incorporating another phase lag caused by a temperature gradient, resulting in the dual-phase-lag (DPL) model. Furthermore, standard electrometric constitutive models are insufficient for describing soft tissue's complex mechanical behavior. A thorough understanding of the thermal [13, 14, 15, 16] and biothermo-mechanical [17, 18, 19] responses of biological tissues will aid in the design and optimization of heat transfer processes in biological tissues. Several models were proposed to study the mechanical behavior of biological tissues [20,21]. [PERSON] and [PERSON] [22] studied the non-linear viscoelastic laws for soft biological tissues. A generic physics-informed neural-network-based constitutive model for soft biological tissues was established by [PERSON] et al. [23]. Furthermore, [PERSON] and [PERSON] [24] studied the non-linear and time-dependent properties of soft biological tissues containing collagen.  
The primary goal of this study was to introduce the hybrid boundary element method (HBEM) scheme for describing nonlinear fractional bio-thermomechanical interactions in functionally graded anisotropic (FGA) soft tissues, which is based on the local radial basis function collocation method (LRBFCM) and GBEM to solve the bioheat governing equation and on the boundary element method (BEM) to solve the poroelastic governing equation. The calculation results were graphed to show the effects of the anisotropic property and graded and fractional parameters on bio-thermal stress. These findings also supported the proposed hybrid scheme's validity and efficiency.