A transcranial photobiomodulation field modeling method based on optical-thermal bidirectional coupling

Abstract

This invention discloses a transcranial photobiological modulation field modeling method based on photothermal bidirectional coupling, applied to the interdisciplinary field of biomedical engineering and applied mathematics. It addresses the problem that current tPBM light field modeling techniques neglect the photothermal bidirectional coupling mechanism and fail to consider the dynamic influence of thermal feedback on the light field distribution, leading to errors in light field prediction and inaccurate safety assessments. This invention establishes a real-time temperature-light absorption coefficient model to accurately characterize the impact of temperature changes on the light field distribution. Furthermore, it combines the Pennes biothermal transfer equation, commonly used to describe biological heat transfer, to model the tissue temperature rise induced by tPBM irradiation, thereby constructing a complete photothermal bidirectional coupling simulation framework. This invention effectively improves the accuracy of light field prediction by dynamically capturing the interaction between the light and temperature fields.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of biomedical engineering and applied mathematics, and specifically relates to a tPBM light field modeling technique. Background Technology

[0002] Transcranial photobiomodulation (tPBM) is an emerging non-invasive neuromodulation technique that uses near-infrared light (typically with wavelengths between 600-1100 nm) to penetrate the skull and irradiate the brain, thereby modulating and repairing neural function. Its core therapeutic mechanism is believed to improve neural function by enhancing brain energy metabolism, showing significant application potential in the treatment of neurodegenerative diseases (such as Alzheimer's and Parkinson's), stroke, and other neurological disorders, as well as in areas such as cognitive enhancement and motor function recovery. However, similar to other non-invasive brain stimulation techniques such as transcranial magnetic stimulation (TMS) and transcranial electrical stimulation (tES), the clinical translation of tPBM faces significant heterogeneity in individual efficacy. This heterogeneity stems primarily from two aspects: first, individual anatomical and biophysical heterogeneity (such as differences in skull thickness, brain sulcus structure, and tissue optical properties among individuals) leads to deviations between the actual and expected stimulation targets; second, tPBM exhibits a unique biphasic dose-response relationship (low energy levels show a promoting effect, while high energy levels may turn into inhibition or even negative effects). Therefore, developing high-fidelity, personalized intracranial light field prediction models has become a key technological bottleneck in achieving precise treatment of tPBM.

[0003] Existing tPBM light field modeling methods mainly fall into two categories, both of which have limitations. Analytical models, driven by computational efficiency, pursue rapid computation by simplifying the geometry and physical processes of biological tissues. For example, Erkol et al. derived analytical solutions to time-dependent diffusion equations based on Bessel functions in simplified 2D circular and 3D cylindrical geometries. While these methods are highly efficient and accurate in regular shapes, they struggle to capture the light propagation process within the complex, layered geometry of a real human head, leading to significant prediction errors. Numerical models, driven by accuracy, overcome the limitations of analytical models by precisely characterizing the complex anatomical geometry of biological tissues to more accurately simulate light propagation. These include Monte Carlo simulations (such as MMC, MCVM, and MCML) and finite element analysis (such as NIRFAST). In recent years, the introduction of technologies such as wide-field illumination source modeling, GPU acceleration, and mesh-based Monte Carlo methods has further improved the computational efficiency and simulation capabilities of numerical models. However, all models are based on a fundamental assumption: the energy of photons absorbed by tissues is either consumed through photochemical reactions or permanently lost from the system. This assumption completely ignores the process of its conversion into heat energy, and the resulting thermal field's influence on the light field distribution, ultimately leading to modeling errors.

[0004] In actual photothermal biomass chemoradiography (tPBM), part of the photon energy absorbed by biological tissue is used for ATP production through biochemical processes (such as absorption by light-sensitive molecules like cytochrome c oxidase, promoting mitochondrial electron transport), while the other part is inevitably converted into heat energy through non-radiative relaxation, leading to a local temperature increase. Experimental measurements show that skin temperature can rise to 40°C under typical tPBM stimulation parameters. The photothermal effect is not a negligible byproduct, but a key factor with both therapeutic and safety significance: it is actively used for tumor ablation in photothermal therapy, and in tPBM, it is a core parameter ensuring treatment safety. Therefore, the forward mathematical physics of photothermal conversion has been extensively studied. Currently, the photothermal conversion effect in biological tissue is mainly described by the Pennes bioheat equation (PBE), and widely used computational models have been developed based on this. Specifically, this energy conversion from light field to temperature field can be viewed as a unidirectionally coupled computational framework: first, optical simulation is performed to calculate the light field distribution and energy deposition; then, the energy deposition is input as a heat source term into the PBE to solve for the temperature field. Although this one-way coupling method can simulate the increase in tissue temperature caused by light energy during tPBM, it does not model the effect of the temperature increase on the light field, and therefore cannot improve the accuracy of light field prediction through this one-way modeling.

[0005] According to the first law of thermodynamics, energy dissipated in the form of heat is not a passive loss, but rather an active alteration of the system's state. Although an increase in temperature does not change the fundamental physical laws governing photon propagation (i.e., the governing equations remain unchanged), it does affect the light field distribution by modulating the optical properties of biological tissues (especially the absorption coefficient). Experimental studies have confirmed that the absorption coefficient of skin tissue exhibits a significant temperature dependence (approximately 10.1% increase for every 1°C increase in temperature). This characteristic is primarily due to the high water content of the tissue and the strong temperature sensitivity of water in the near-infrared spectral region. Sfareni et al. further discovered that when the temperature rises from 20°C to 40°C, the water absorption coefficient increases by approximately 10%, and the absorption peak exhibits a blue shift (from 971 nm to 966 nm). Within the physiological temperature range (≤42°C), the absorption coefficients of human skin and brain tissue show an approximately linear positive correlation with temperature. This correlation between temperature and absorption coefficient provides a theoretical basis for constructing a reverse coupling mechanism between the temperature field and the light field.

[0006] In summary, the core problem with current tPBM optical field modeling technology is that it ignores the photothermal bidirectional coupling mechanism and does not consider the dynamic impact of thermal feedback on the optical field distribution, resulting in optical field prediction errors and inaccurate safety assessments. Summary of the Invention

[0007] To address the aforementioned technical problems, this invention proposes a transcranial photobiological modulation field modeling method based on photothermal bidirectional coupling. By dynamically capturing the interaction between the light field and the temperature field (i.e., the bidirectional feedback process in which the light field influences the temperature field and the temperature field modulates the light field in turn), the accuracy of light field prediction is improved.

[0008] The technical solution adopted in this invention is: a transcranial photobiological modulation field modeling method based on photothermal bidirectional coupling, comprising:

[0009] S1. Construct a light field distribution model within the head tissue, including mathematical description, weakening, discretization, and linear system construction of the light transmission process;

[0010] S2. Based on the light field distribution model constructed in step S1, construct a temperature field distribution model within the head tissue, including the mathematical and physical description, weakening, discretization, and linear system construction of the heat conduction process;

[0011] S3. Construct a real-time temperature-dependent light absorption coefficient model to introduce dynamic feedback of temperature on optical parameters into the light field distribution model constructed in step S1.

[0012] S4. Construct the head model and solve the linear system of the temperature field distribution model under the head model to obtain the temperature, luminous flux and light absorption coefficient of the current iterative light-irradiated area.

[0013] The beneficial effects of this invention: Existing tPBM light field simulation algorithms fail to consider the feedback effect of temperature on the light field, leading to decoupling errors between the simulated light field and temperature field. Specifically, this manifests as an overestimation of intracranial light flux and an underestimation of tissue temperature. To address this bias, this invention proposes a transcranial light-modulated field modeling method based on photothermal bidirectional coupling. By establishing a real-time temperature-light absorption coefficient model (RT-ACM), the influence of temperature changes on the light field distribution is accurately characterized. Furthermore, by combining the Pennes biological heat transfer equation (PBE), commonly used to describe biological heat transfer, the tissue temperature rise induced by tPBM irradiation is modeled, thus constructing a complete photothermal bidirectional coupling simulation framework. This method, on the one hand, overcomes the physical simplification assumptions in traditional single-field modeling, improving the physical realism of the light field simulation; on the other hand, by introducing dynamic feedback of temperature on optical parameters, it effectively compensates for the decoupling errors caused by neglecting photothermal interactions in traditional PBE temperature modeling. Compared with the existing uncoupled modeling methods mentioned in the background art, the coupling strategy proposed in this invention is expected to significantly improve the simulation accuracy of the light field and temperature field in the tPBM process, and provide a more reliable theoretical basis and numerical support for the formulation of personalized and precise tPBM stimulation schemes. Attached Figure Description

[0014] Figure 1 This is a flowchart of a transcranial photobiological modulation field modeling method based on photothermal coupling according to the present invention.

[0015] Figure 2 This is a schematic diagram illustrating the verification results of the analytical solution based on the regular cylindrical model in an embodiment of the present invention. Detailed Implementation

[0016] This invention establishes the laws governing light propagation and heat diffusion in biological tissues, as well as the forward coupling of the light field to the temperature field, based on the Photon Diffusion Equation (PDE) and the Photon Behavior Equation (PBE), respectively. Then, by constructing a real-time temperature-dependent absorption coefficient model (RT-ACM), the influence of temperature changes on the absorption coefficient of head tissue is characterized, and the reverse coupling of the thermal field to the light field is constructed, forming a photo-thermal bidirectional coupling closed-loop framework (PT-BCM).

[0017] The method of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0018] like Figure 1The flowchart shown is a modeling method for transcranial photobiological modulation field based on photothermal bidirectional coupling according to the present invention. The specific steps are as follows:

[0019] S1. Construct a light field distribution model within the head tissue, including mathematical description, weakening, discretization, and linear system construction of the light transmission process.

[0020] (1) Mathematical and physical description of the optical transmission process

[0021] When the isotropic fluence significantly exceeds the directional flux, the optical field is considered "diffuse." This occurs in regions dominated by scattering, far from the source and boundaries, provided that the fluence does not change rapidly over time (e.g., on a sub-picosecond scale). This assumption allows for simplification from the general radiative transfer equations describing anisotropic fields to a diffusion approximation applicable to isotropic fluence.

[0022]

[0023] In the formula, The solution domain for the light diffusion equation is represented as the head tissue space. Represents the spatial gradient; Indicates the imaginary part; Indicates the location at coordinate point The control frequency is Luminous flux, measured in units of: ; Indicates the diffusion coefficient; and They respectively represent the locations The absorption coefficient and the reduced scattering coefficient; Indicates that it is located at the speed of light; This represents the isotropic source term (or the distribution of light sources).

[0024] A Robin boundary condition with refractive index mismatch was applied to the outer boundary. This condition stipulates that the outgoing flux equals the local flux multiplied by a coefficient simulating internal reflection, thus preventing light from flowing back from the air into the tissue. This relationship is defined by the following equation:

[0025]

[0026] In the formula, Indicates the luminous flux on the scalp surface; Indicates the surface of the scalp; Represents the coordinates of points on the scalp surface; This represents the normal vector pointing outwards from the scalp surface; The degree of mismatch in refractive index between the scalp and air can be obtained using Fresnel's theorem, as follows:

[0027]

[0028] In the formula, Indicates that photons travel from a refractive index of 1000 to 10000. area Entering the refractive index The critical angle at which total internal reflection occurs when exposed to external air; This represents the vertical incident reflectivity at the boundary. Typically, the refractive index of the external air is set to the vacuum refractive index, i.e. .

[0029] (2) Weakening

[0030] Based on the mathematical description of the light propagation process, we will now use the finite element method to discretize and approximate the solution of the complex geometric space. The first step of the finite element method is weakening, that is, using the test function of the first-order Sobolev space. Weakening formula (1).

[0031]

[0032] in, Represents a first-order Sobolev space;

[0033] Further apply the boundary conditions shown in formula (2) to constrain the process, i.e., replace... Part of, by

[0034]

[0035] (3) Discretization

[0036] The weakened formula (5) is further spatially discretized. The discretization method depends on the type of discrete mesh (tetrahedral or hexahedral) of the solution domain; different mesh types correspond to different shape functions. Given mesh parameters... , Represents a discrete region. The number of grid nodes. Approximate discrete photon flux. It can be represented as:

[0037]

[0038] In the formula, Represents discrete nodes Luminous flux at the location; It is a three-dimensional linear shape function.

[0039] Substituting the weak scheme formula (4) of the optical transmission mathematical equation into the above discretized space We can obtain:

[0040]

[0041] By further changing the order of integration and summation, we obtain the finite element form of the following mathematical equations for optical transmission.

[0042]

[0043] (4) Construction of linear systems

[0044] Rearranging the above integral equation (8), we obtain the following system of linear equations:

[0045]

[0046] In the formula, These are terms related to geometric and tissue physical properties. Represents a real number of N*N dimensions; It is an item related to the tPBM irradiation source; It is a numerical solution obtained from the system of equations of this linear system, which serves as an ideal approximation of the photon flux. , , and The calculation is as follows:

[0047]

[0048]

[0049]

[0050]

[0051] S2. Using the light field distribution constructed in S1, construct a temperature field distribution model, including the mathematical and physical description of the heat conduction process, weakening, discretization, and linear system construction.

[0052] (1) Mathematical and physical description of heat conduction

[0053] The energy conservation of biological heat transfer is described by the biological heat conduction equation, which is derived from Fourier's law of heat conduction. Therefore, this equation can be used to directly calculate the temperature changes caused by the absorption of light in different layers of biological tissue.

[0054]

[0055] In the formula, Indicates the temperature of the head tissue; , and These represent tissue density, tissue specific heat capacity, and tissue thermal conductivity, respectively. Indicates the core body temperature; and This indicates blood perfusion rate and specific heat capacity. This indicates the heat energy contribution from blood perfusion; The unit for expressing the thermal energy contribution from light energy is W / m³.

[0056] Taking into account the heat exchange that exists between the scalp surface and the air, the Robin boundary is used to describe this heat exchange on the scalp surface.

[0057]

[0058] In the formula, The heat convection coefficient is set as follows: ; The ambient temperature is set to 25 °C.

[0059] (2) Weakening

[0060] The weakening of formula (13) is achieved by multiplying by the test function. And then further perform distribution integration to obtain the result.

[0061]

[0062] Further apply the boundary conditions shown in Equation (15), i.e., replace In part, a weak scheme of the mathematical physics equation for heat conduction is obtained.

[0063]

[0064] (3) Discretization

[0065] For discrete space Discrete approximation of temperature This can be expressed as:

[0066]

[0067] In the formula, Represents discrete nodes Temperature at that location; Represents a three-dimensional linear shape function.

[0068] Substituting the weak scheme of the heat conduction mathematical equation (Equation (17)), we obtain the following discretized form:

[0069]

[0070] By changing the order of the integral and summation operators, we obtain the following finite element scheme for PBE.

[0071]

[0072] (4) Construction of linear systems

[0073] Reorganizing the discrete form of the above mathematical physics equations for heat conduction into a system of linear equations, we obtain the following ordinary differential linear system:

[0074]

[0075] In the formula, These are the mass matrix and the stiffness matrix, respectively; Let each of the two vectors be related to the source term, where... Related to constant heat exchange caused by ambient temperature and blood perfusion, Related to the thermal conversion caused by the light source; This represents the discrete temperature solution to be found for the above system. , , and The calculation is as follows:

[0076]

[0077]

[0078]

[0079]

[0080] S3. Construct a temperature feedback model and complete the construction of the photothermal bidirectional coupling framework.

[0081] To establish the reverse coupling relationship, we further developed a real-time temperature-dependent absorption coefficient model (RAM). Specifically, based on existing experimental studies, we determined the temperature sensitivity coefficients of the scalp, skull, and brain at a wavelength of 810 nm. The coefficients for each tissue are assigned as follows: (1) Scalp ( (2) The coefficient is calculated by formula (25); (3) The coefficient of the skull is set to zero because of its low water content and extremely small thermal response; (4) The coefficient of the brain tissue is 1.865×10⁻ from previous studies. 6 (mm·℃)⁻¹ is used as the temperature sensitivity coefficient.

[0082]

[0083] In the formula, , This represents the temperature sensitivity of the absorption coefficient of different chromophores, where the temperature sensitivity of the absorption coefficient of the HHb chromophore is... Set as , Thermosensitive absorption coefficient of chromophore Set as ;water( The absorption coefficient of chromophores is temperature sensitive. Set as ; and These represent oxygen saturation and water content, respectively.

[0084] The absorption coefficients of each tissue in the model are updated by substituting the aforementioned RT-ACM model into the light field distribution model constructed in step S1. The optical-thermal bidirectional coupling framework was completed.

[0085] S4, bidirectional coupled framework solution and time update, including head model construction and solution of ordinary differential equations.

[0086] (1) Head model construction

[0087] The head model construction consists of two parts: segmentation and mesh construction. First, based on the subject's MRI (Magnetic Resonance Imaging), the head domain is divided into different tissue domains, which should at least include three layers of brain tissue: scalp, skull, and brain. Then, the segmented tissue domains are discretized, usually using tetrahedral or hexahedral meshes. Different discretized meshes require different shape functions for description.

[0088] (2) Solving the system of ordinary differential equations

[0089] The ordinary differential equation system (21) is solved discretized in the time domain. Discretization methods include the conditionally stable Euler method and Runge-Kutta method, as well as the unconditionally stable Crank-Nicolson method. Because the Crank-Nicolson method is unconditionally stable, it can effectively reduce the influence of mesh size on numerical computation. The Crank-Nicolson method is used as an example for demonstration below, where the mesh size is typically... .

[0090]

[0091] In the formula, Indicates the current iteration step. This indicates the next iteration step.

[0092] In this embodiment, step S5 is also included, which is used to evaluate the model error.

[0093] To verify the proposed photothermal bidirectional coupling model, we compared its numerical solution with the numerical and analytical solutions of the traditional uncoupled model, and used the relative difference measure (RDM) to quantify the difference, which is the model error between the numerical and analytical solutions.

[0094]

[0095] In the formula, This represents the Euclidean norm. The temperature, luminous flux, and absorption coefficient in the irradiated area obtained by using the method of the present invention as a control group are represented by a column vector of the number of grid nodes covered by the irradiated area (the cylindrical area below the irradiated area). This indicates that the existing uncoupled model in the background technology is used as the reference group of column vectors.

[0096] Figure 2 This is a schematic diagram illustrating the verification results of the analytical solution based on the regular cylindrical model in this embodiment of the invention. The diagram shows a benchmark comparison between the analytical solution and the numerical results of the coupled / uncoupled models, including specific settings and comparison results: Sub-figure (a) presents the computational domain and mesh generation used in the numerical simulation; sub-figure (b) compares the temperature distribution on the inner cylinder surface using the three methods and gives the RDM values ​​of each numerical model relative to the analytical reference solution. The consistent range of color bars between the analytical and numerical solutions initially indicates that the calculation results have physical rationality and numerical consistency. Quantitatively, the RDM values ​​of both models decrease exponentially with irradiation time, which is consistent with the inherent error diffusion characteristics of the diffusion equation, proving the mathematical consistency of the numerical solver used in this invention. Compared to the uncoupled model, the coupled model consistently maintains a lower RDM value, with a relative error reduction of 2.35%–8.62%, thus establishing its superior prediction accuracy throughout the simulation period. Simultaneously, it can be seen that the uncoupled model exhibits a larger spatial deviation, particularly noticeable near the geometric boundaries, indicating that the numerical error of the uncoupled model mainly occurs in shallow tissues. Figure 2 In the table, Height represents the height, Inner represents the inner diameter, Outer represents the outer diameter, Sources represents the light source settings for the analytical solution, Analytic represents the numerical results of the analytical model, Coupled represents the numerical results of the coupled model, Uncoupled represents the numerical results of the uncoupled model, and Time represents the time.

[0097] In summary, this study developed a novel photothermal bidirectional coupling model to simulate the spatiotemporal dynamics of photo-thermal coupling in transcranial photobiological modulation. This model overcomes the fundamental limitations of traditional uncoupled methods by combining the Pennes biothermal equation with a real-time temperature-light absorption coefficient model, forming a closed photothermal feedback loop. The proposed method ensures model consistency with the law of energy conservation, enabling predictions that provide a physically consistent explanation for attenuated photon flux and the resulting enhanced local energy deposition. Therefore, the model has the potential to correct the inherent systematic overestimation of photon flux and underestimation of tissue temperature errors in uncoupled methods. Numerical validation confirmed the accuracy of the coupled equation solution. Furthermore, this framework enables spatiotemporal dynamic simulation of light field modeling through photothermal coupling, providing a feasible tool for studying tPBM modulation modes with contrasting temporal modes (such as continuous wave modulation and pulsed wave modulation). In conclusion, this photothermal bidirectional coupling model provides a computational framework integrating photothermal physical processes, laying the foundation for more accurate dose calculations and becoming an important tool for mechanistic research.

[0098] Those skilled in the art will recognize that the embodiments described herein are for the purpose of helping to understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Various modifications and variations can be made to the invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the invention should be included within the scope of the claims of the invention.

Claims

1. A transcranial photobiological modulation field modeling method based on photothermal bidirectional coupling, characterized in that, include: S1. Construct a light field distribution model within the head tissue, including constructing a mathematical description of the light transmission process and discretizing the space of complex geometry based on the finite element method. The solution process includes: weakening and discretizing the mathematical description of the light transmission process, and constructing a linear system based on the discretization results. S2. Based on the light field distribution model constructed in step S1, construct a temperature field distribution model within the head tissue, including constructing a mathematical and physical description of the heat conduction process, and discretizing and solving the complex geometric space based on the finite element method. The solution process includes: weakening and discretizing the mathematical and physical description of the heat conduction process, and constructing a linear system based on the discretization results. S3. Construct a real-time temperature-dependent light absorption coefficient model to introduce dynamic feedback of temperature on optical parameters into the light field distribution model constructed in step S1. S4. Construct the head model and solve the linear system of the temperature field distribution model under the head model to obtain the temperature, luminous flux and light absorption coefficient of the current iterative light-irradiated area.

2. The transcranial photobiological modulation field modeling method based on photothermal bidirectional coupling according to claim 1, characterized in that, Step S1 involves constructing the mathematical description of the optical transmission process as follows: Assuming that the flux does not change rapidly over time, the optical field is considered "diffuse" when the isotropic flux significantly exceeds the directional flux; this occurs in regions dominated by scattering, far from the source and boundaries. Based on this assumption, the general radiative transfer equation describing anisotropic fields simplifies to a diffusion approximation applicable to isotropic flux: ； in, The solution domain for the light diffusion equation is represented as the head tissue space. Represents the spatial gradient; Indicates the imaginary part; Represents the coordinates of the points located on the scalp surface. The control frequency is luminous flux; Indicates the diffusion coefficient. ; and These represent the coordinates of the points located on the scalp surface. The absorption coefficient and the reduced scattering coefficient; Represents the coordinates of the points located on the scalp surface. the speed of light; Represents the coordinates of the points located on the scalp surface. The control frequency is isotropic source terms; Apply a refractive index mismatch Robin boundary condition to the outer boundary, which stipulates that the outgoing flux is equal to the local flux multiplied by a coefficient simulating internal reflection.

3. The transcranial photobiological modulation field modeling method based on photothermal bidirectional coupling according to claim 2, characterized in that, Step S2 involves constructing a mathematical and physical description of the heat conduction process as follows: The energy conservation of biological heat transfer is described by the biological heat conduction equation, which is derived from Fourier's law of heat conduction. Therefore, this equation can be used to directly calculate the temperature changes caused by the absorption of light in different layers of biological tissue. ； in, Indicates the temperature of the head tissue; , and These represent tissue density, tissue specific heat capacity, and tissue thermal conductivity, respectively. Indicates the core body temperature; and This indicates blood perfusion rate and blood specific heat capacity; This indicates the heat energy contribution from blood perfusion; The unit representing the thermal energy contribution from light energy is: Considering the heat exchange between the scalp surface and the air, the heat exchange on the scalp surface is described using a natural convection boundary.

4. The transcranial photobiological modulation field modeling method based on photothermal bidirectional coupling according to claim 3, characterized in that, The real-time temperature-dependent light absorption coefficient model constructed in step S3 is as follows: The temperature sensitivity coefficients of the scalp, skull, and brain tissue at selected wavelengths are determined. The formula for calculating the temperature sensitivity coefficient of the scalp is as follows: ； in, Indicates different chromophores The absorption coefficient is temperature sensitive. , and These represent oxygen saturation and water content, respectively.

5. The transcranial photobiological modulation field modeling method based on photothermal bidirectional coupling according to claim 4, characterized in that, Step S3 introduces dynamic feedback of temperature on optical parameters into the light field distribution model constructed in step S1. Specifically, it replaces the absorption coefficients of the scalp, skull, and brain tissue at selected wavelengths with the temperature sensitivity coefficients of the scalp, skull, and brain tissue in the light field distribution model constructed in step S1.

6. The transcranial photobiological modulation field modeling method based on photothermal bidirectional coupling according to claim 5, characterized in that, The head model construction described in step S4 is as follows: the head domain is divided into different tissue domains based on the subject's MRI, including at least three layers of brain tissue: scalp, skull and brain; then, the segmented tissue domains are discretized.

7. The transcranial photobiological modulation field modeling method based on photothermal bidirectional coupling according to claim 6, characterized in that, Discretization is performed using tetrahedral or hexahedral meshes.

8. The transcranial photobiological modulation field modeling method based on photothermal bidirectional coupling according to claim 7, characterized in that, In step S4, the linear system of the temperature field distribution model under this head model is solved by using one of the following methods: the conditionally stable Euler method, the conditionally stable Runge-Kutta method, or the unconditionally stable Crank-Nicolson method.

Images