Method and system for constructing a patient full parameter digital twin and performing treatment simulation

Abstract

This invention belongs to the field of computational medicine, specifically relating to a method and system for constructing a fully parameterized digital twin of a patient and simulating treatment. It aims to address the problems of existing digital twin models, such as their inability to integrate multi-scale data, lack of dynamic optimization capabilities, and disconnect from clinical decision-making processes. The invention includes: acquiring and integrating multi-source heterogeneous patient data; performing cross-scale fusion modeling based on the data to construct a personalized initial parameter set for the digital twin, including genetic background correction; driving a multiphysics coupling engine to simulate intervention effects and constructing an inverse optimization problem with real-time monitoring data as dynamic constraints to solve for the optimal intervention parameters; performing a validation simulation of the optimized plan and generating a treatment report including quantitative assessment and risk warning. By constructing a fully parameterized, mechanistic model and introducing a dynamic closed-loop optimization validation process, this invention significantly improves the personalized accuracy, safety, and clinical decision support value of treatment plan simulation.

Description

Technical Field

[0001] This invention belongs to the field of computational medicine, specifically relating to a method and system for constructing a full-parameter digital twin of a patient and simulating treatment. Background Technology

[0002] Digital twin technology, serving as a crucial bridge connecting the physical world and virtual space, has been widely applied in the industrial sector. In recent years, its application in constructing digital twins of patients for personalized treatment planning and outcome prediction has become a cutting-edge direction in precision medicine. Existing related technologies mainly suffer from the following limitations:

[0003] First, there are shortcomings in the comprehensiveness and depth of model construction. Most existing solutions rely on single imaging data (such as CT and MRI) to construct anatomical models, or combine some time-series physiological signals (such as ECG and blood pressure) for simple dynamic fitting. These models fail to systematically integrate data from macroscopic anatomy and systems physiology to microscopic molecular and genetic data, resulting in digital twins that are only superficial or localized models, unable to reflect the complete biological mechanisms behind diseases. In particular, the lack of characterization of individual physiological responses based on genetic background fundamentally limits the personalization and predictive accuracy of the models.

[0004] Secondly, in terms of treatment simulation and optimization, existing methods often remain at the open-loop simulation level. That is, they simulate and predict the effect of a single treatment based on fixed model parameters. However, in real clinical scenarios, the patient's physiological state is dynamically changing, and the optimal treatment plan needs to be dynamically adjusted based on the patient's real-time feedback. Existing technologies lack a closed-loop mechanism that can use real-time monitoring data during the treatment process as dynamic constraints to optimize treatment parameters in reverse. This makes it difficult to directly apply the simulated ideal treatment plan to dynamically changing real patients, and it cannot achieve real-time personalized correction of the treatment plan.

[0005] Furthermore, existing technical solutions suffer from loose coupling between different stages, lacking a logical closed loop. Anatomical modeling, physiological parameterization, treatment simulation, and treatment optimization are typically treated as independent modules. Barriers exist in the transfer of data and parameters across different scales, failing to form a complete technical closed loop from data to model, from simulation to optimization, and from optimization feedback to validation. This results in the final treatment recommendations lacking rigorous validation based on the same digital twin model, significantly reducing their reliability and clinical guidance value.

[0006] Therefore, the urgent technical problem to be solved is: how to construct a full-parameter digital twin system for patients that can integrate multi-scale biomedical features, possess high-fidelity physiological mechanisms, and introduce real-time dynamic constraints for closed-loop optimization and verification, so as to achieve truly personalized, dynamic, and verifiable treatment plan simulation and deduction. Summary of the Invention

[0007] To address the aforementioned problems in existing technologies, namely, the inability of existing digital twin models to integrate multi-scale data, lack of dynamic optimization capabilities, and disconnect from clinical decision-making processes, this invention provides a method for constructing a full-parameter digital twin of a patient and simulating treatment, the method comprising:

[0008] Acquire and integrate static anatomical data, dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic characteristic data of target patients to construct a multi-source heterogeneous biomedical feature pool.

[0009] The static anatomical structure data is analyzed to generate a three-dimensional anatomical structure model, and combined with the dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic characteristic data to perform cross-scale fusion modeling, forming the initial parameter set of the personalized digital twin;

[0010] The personalized digital twin initial parameter set drives the multi-physics field coupled digital twin kernel engine to simulate the effect of preset intervention measures on the patient's full-parameter digital twin and output a set of predicted physiological state parameters.

[0011] An inverse optimization problem is constructed based on the difference between the preset clinical goals and the predicted physiological state parameter set. Dynamic physiological state constraints based on real-time patient monitoring data are introduced during the optimization process to solve and generate optimized intervention plan parameters.

[0012] The optimized intervention parameters are input into the digital twin kernel engine for verification simulation, and a treatment simulation report is generated based on the set of verification physiological state parameters output by the simulation.

[0013] Furthermore, the static anatomical structure data is analyzed to generate a three-dimensional anatomical structure model, the method of which is as follows:

[0014] Spatial registration and voxel-level semantic fusion of multimodal static anatomical imaging data are performed to generate unified 3D volume data with tissue category annotations.

[0015] The three-dimensional volume data is subjected to hierarchical segmentation based on prior anatomical atlases, and the boundary surfaces of different anatomical structures and substructures are extracted and reconstructed to form a three-dimensional mesh model that expresses the geometric morphology.

[0016] Based on the aforementioned three-dimensional mesh model, a multi-scale topological map describing the spatial adjacency relationships between organs and the physiological connectivity within tissues is analyzed and constructed.

[0017] The material property parameters characterizing the biophysical properties of the tissue are assigned to the corresponding units of the three-dimensional mesh model according to the tissue category label and the multi-scale topology map, thereby generating a computable three-dimensional anatomical structure model that integrates geometric, topological and physical properties. The three-dimensional anatomical structure model provides the digital twin kernel engine with a discretized spatial domain for multi-physics coupling calculation.

[0018] Furthermore, by combining the dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic characteristic data, cross-scale fusion modeling is performed to form the initial parameter set for the personalized digital twin. The method is as follows:

[0019] Based on the dynamic physiological function data and the time-series metabolic parameter data, a system state matrix is ​​constructed, and a system identification algorithm is applied to solve the set of ordinary differential equations describing macroscopic physiological behavior to obtain a primary set of kinetic parameters.

[0020] The molecular and genetic feature data are mapped to a pre-set biological pathway knowledge graph. By calculating the centrality measure of specific biomolecular entities in the graph topology and the weights of their associated pathways, a set of genetic regulatory coefficients are generated.

[0021] Using the genetic regulation coefficient, the primary dynamic parameter set is corrected through a preset parameter modulation function. The modulation function takes the genetic regulation coefficient as input and performs element-wise transformation on the state transition matrix of the primary dynamic parameter set to obtain the corrected system parameters under genetic background constraints.

[0022] A cross-scale mapping model is established with physical field smoothness and conservation law as joint optimization objectives. The modified system parameters are used as the source field and mapped to the physical attribute target field of each discrete unit of the three-dimensional anatomical structure model to generate the initial parameter set of the personalized digital twin.

[0023] Furthermore, the personalized digital twin initial parameter set drives a multi-physics coupled digital twin kernel engine to simulate the effect of preset intervention measures on the patient's full-parameter digital twin, outputting a predicted set of physiological state parameters. The method is as follows:

[0024] Based on the mesh topology of the three-dimensional anatomical structure model, the blood rheological parameters and tissue mechanical modulus in the initial parameter set of the personalized digital twin are mapped to fluid domain units and solid domain nodes, respectively, to construct a multiphysics calculation model containing non-uniform material properties.

[0025] Based on the multiphysics calculation model, a simultaneous solution mechanism including the Navier-Stokes equations, nonlinear constitutive equations, and reaction-convection-diffusion equations is established. The preset intervention measures are quantified as boundary conditions in the solution mechanism. For drug treatment, an inlet concentration function or mass source term is defined, and for physical intervention, a contact boundary or displacement constraint is defined.

[0026] The simultaneous solution mechanism is driven to perform fluid-structure interaction iteration within the time dispersive walk, transferring the shear force and pressure on the fluid side to the solid side to calculate the deformation, and feeding back the deformation displacement to the fluid side to update the mesh, until the full spatiotemporal numerical simulation is completed;

[0027] For the physical field evolution data obtained from the simulation, a virtual sampling probe is set in the anatomical coordinate system, and the time series and statistical features of key variables are extracted by shape function interpolation, and the predicted physiological state parameter set is compiled and generated.

[0028] Furthermore, the multiphysics computational model is constructed as follows:

[0029] The three-dimensional anatomical model is segmented based on physical properties to identify and mark the vascular lumen region, vascular wall region, and surrounding target tissue region.

[0030] Based on the results of the region segmentation, tetrahedral or hexahedral mesh units in the three-dimensional anatomical structure model are assigned region labels, and based on the parameter subset associated with each region label in the personalized digital twin initial parameter set, specific material property values ​​are assigned to each mesh unit or node, wherein the property values ​​of the vascular lumen unit are set according to the blood rheological parameters, and the property values ​​of the vascular wall and target tissue units are set according to the tissue mechanical modulus.

[0031] For mesh cells or nodes located at the boundaries of different material regions, a distance-based linear or nonlinear interpolation algorithm is used to smooth the material property values ​​and generate the multiphysics calculation model containing non-uniform material properties.

[0032] Furthermore, the simultaneous solution mechanism is driven to perform bidirectional fluid-structure interaction iteration within the time dispersive walk, and the method is as follows:

[0033] A partitioned, strongly coupled iterative strategy is adopted to alternately solve the fluid dynamics and solid mechanics control equations and exchange boundary data in each global time step.

[0034] With the current solid domain mesh shape fixed, solve the Navier-Stokes equations to obtain the pressure and velocity fields of the fluid domain, and calculate the fluid loads at the fluid-solid interface;

[0035] The fluid load is applied to the corresponding boundary node of the solid domain, and the deformation displacement field of the solid domain is obtained by solving the nonlinear constitutive equation.

[0036] Based on the deformation displacement field, the node coordinates of the fluid domain mesh are updated using an arbitrary Lagrange-Euler method;

[0037] Determine whether the forces and displacements at the fluid-structure interface meet the preset convergence conditions. If not, repeat the alternating solution based on the updated mesh and load until convergence, and then proceed to the next time step.

[0038] Furthermore, an inverse optimization problem is constructed based on the difference between the preset clinical goals and the predicted physiological state parameter set. Dynamic physiological state constraints based on real-time patient monitoring data are introduced during the optimization process to solve for and generate optimized intervention plan parameters. The method is as follows:

[0039] Based on the quantitative difference between the predicted physiological state parameter set and the preset clinical target, an inverse optimization problem is constructed with minimizing the difference as the objective function;

[0040] The real-time patient monitoring data is converted into a feasible domain boundary condition that evolves over time, and this boundary condition is introduced as an inequality constraint term into the inverse optimization problem to form a dynamic constraint optimization model.

[0041] The dynamic constraint optimization model is solved iteratively using a gradient-based optimization algorithm. In each iteration, the boundary conditions of the feasible region are updated synchronously to reflect the latest monitoring data. The search direction is guided by calculating the constraint violation penalty term until the algorithm converges.

[0042] The decision variable values ​​corresponding to the convergence are output as parameters of the optimized intervention scheme.

[0043] Furthermore, the real-time patient monitoring data is transformed into a feasible domain boundary condition that evolves over time, and this boundary condition is introduced as an inequality constraint term into the inverse optimization problem. The method is as follows:

[0044] The patient's real-time monitoring data is processed in real time, and the trend characteristics and fluctuation range of key physiological parameters are extracted through a sliding time window;

[0045] Based on the extracted trend features and fluctuation range, a dynamic feasible region describing the safety boundary of physiological state is established, and this dynamic feasible region is modeled as a time-varying inequality constraint of decision variables in the inverse optimization problem.

[0046] During the iterative solution process of the optimization algorithm, the specific parameters of the time-varying inequality constraint are updated according to the latest window of the monitoring data to ensure that the parameters of the generated optimized intervention plan always meet the patient's current physiological state boundary.

[0047] Furthermore, the optimized intervention parameters are input into the digital twin kernel engine for verification simulation, and a treatment simulation report is generated based on the set of verification physiological state parameters output by the simulation. The method is as follows:

[0048] The optimized intervention scheme parameters are configured as the input boundary conditions or source terms of the digital twin kernel engine, driving the kernel engine to perform a verification numerical simulation based on the same physical mechanism as the original simulation but with optimized parameters, and outputting a set of verification physiological state parameters containing complete spatiotemporal evolution data.

[0049] The set of verified physiological state parameters is subjected to post-processing analysis, including calculating its achievement degree with the preset clinical goal, the improvement magnitude with the predicted set of physiological state parameters, and assessing its long-term evolution trend and potential risks based on the kinetic stability criterion.

[0050] The results of the post-processing analysis, together with the parameters of the optimized intervention plan and the visualization data of the key simulation process, are structurally integrated according to a preset template to generate the treatment simulation report containing quantitative assessment, plan comparison and risk warning information.

[0051] In a second aspect, the present invention provides a system for constructing a full-parameter digital twin of a patient and performing treatment simulation, for executing a method for constructing a full-parameter digital twin of a patient and performing treatment simulation, the system comprising:

[0052] The data integration module is configured to acquire and integrate static anatomical structure data, dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic characteristic data of the target patient to construct a multi-source heterogeneous biomedical feature pool.

[0053] The digital twin modeling module is configured to parse the static anatomical structure data to generate a three-dimensional anatomical structure model, and combine the dynamic physiological function data, time-series metabolic parameter data and molecular and genetic feature data to perform cross-scale fusion modeling to form the personalized digital twin initial parameter set.

[0054] The multi-physics coupling simulation module is configured to drive the multi-physics coupling digital twin kernel engine with the personalized digital twin initial parameter set, simulate the effect of preset intervention measures on the patient's full-parameter digital twin, and output a set of predicted physiological state parameters.

[0055] The dynamic constraint optimization module is configured to construct an inverse optimization problem based on the difference between the preset clinical goal and the predicted physiological state parameter set, and introduce dynamic physiological state constraints based on real-time monitoring data of patients during the optimization process to solve and generate optimized intervention plan parameters.

[0056] The verification and report generation module is configured to input the parameters of the optimized intervention plan into the digital twin kernel engine for verification simulation, and generate a treatment simulation report based on the set of verification physiological state parameters output by the simulation.

[0057] The beneficial effects of this invention are:

[0058] This invention systematically integrates data from macroscopic anatomy and systems physiology to microscopic molecular and genetic data, and utilizes genetic information to specifically modify physiological dynamics models. The resulting digital twin model overcomes the limitations of traditional methods that rely on population average parameters. This model can intrinsically reflect the anatomical structure and physiological and pathological mechanisms of a specific patient, thereby making more realistic mechanistic predictions of responses to drug treatments, physical interventions, and other interventions, providing a high-fidelity virtual experimental platform for precision medicine.

[0059] Unlike static planning, this invention innovatively introduces dynamic physiological state constraints defined by real-time monitoring data streams during the treatment plan optimization process. This ensures that the intervention parameters (such as drug dosage and surgical parameters) generated by the optimization algorithm are always within the safe boundaries of the patient's current physiological tolerance, enabling the treatment plan to adaptively adjust according to fluctuations in the patient's condition, greatly enhancing the safety and clinical feasibility of the treatment.

[0060] This invention performs a validation simulation of the optimized scheme within the same high-fidelity simulation engine, assesses long-term effects and potential risks, and ultimately generates a structured decision report. This closed-loop process overcomes the limitations of traditional open-loop prediction, providing clinicians with comprehensive decision support that includes quantitative effect comparisons, risk assessments, and visualized evidence. This ensures that the output of digital twins is no longer isolated predictive data, but rather a reliable clinical action reference that has undergone internal cross-validation. Attached Figure Description

[0061] Other features, objects, and advantages of this application will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings:

[0062] Figure 1 This is a flowchart illustrating a method for constructing a full-parameter digital twin of a patient and simulating treatment according to the present invention.

[0063] Figure 2 This is a schematic diagram of the connection relationship of a system for constructing a full-parameter digital twin of a patient and simulating treatment according to the present invention. Detailed Implementation

[0064] The present application will now be described in further detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are for illustrative purposes only and are not intended to limit the invention. Furthermore, it should be noted that, for ease of description, only the parts relevant to the invention are shown in the accompanying drawings.

[0065] It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other. This application will now be described in detail with reference to the accompanying drawings and embodiments.

[0066] The first embodiment of the present invention provides a method for constructing a full-parameter digital twin of a patient and performing treatment simulation, the method comprising:

[0067] Step S10: Acquire and integrate the static anatomical structure data, dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic characteristic data of the target patient to construct a multi-source heterogeneous biomedical characteristic pool;

[0068] Step S20: The static anatomical structure data is parsed to generate a three-dimensional anatomical structure model, and cross-scale fusion modeling is performed by combining the dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic characteristic data to form the initial parameter set of the personalized digital twin.

[0069] Step S30: Drive the multi-physics field coupled digital twin kernel engine with the personalized digital twin initial parameter set to simulate the effect of preset intervention measures on the patient's full-parameter digital twin and output a set of predicted physiological state parameters.

[0070] Step S40: Construct an inverse optimization problem based on the difference between the preset clinical goal and the predicted physiological state parameter set, and introduce dynamic physiological state constraints based on real-time monitoring data of the patient during the optimization process to solve and generate optimized intervention plan parameters.

[0071] Step S50: Input the optimized intervention plan parameters into the digital twin kernel engine for verification simulation, and generate a treatment simulation report based on the set of verification physiological state parameters output by the simulation.

[0072] To more clearly illustrate the method of constructing a full-parameter digital twin of a patient and performing treatment simulation according to the present invention, the following is in conjunction with... Figure 1 The steps in the embodiments of the present invention are described in detail below:

[0073] Step S10: Acquire and integrate the static anatomical structure data, dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic characteristic data of the target patient to construct a multi-source heterogeneous biomedical characteristic pool;

[0074] Specifically, step S10, which involves acquiring and integrating the static anatomical structure data, dynamic physiological function data, temporal metabolic parameter data, and molecular and genetic characteristic data of the target patient to construct a multi-source heterogeneous biomedical feature pool, is executed as follows: First, in acquiring static anatomical structure data, the target patient's medical image files are read using the standard digital medical image transmission protocol, DICOM 3.0. This includes computed tomography images with a slice thickness of no more than 0.625 mm, as well as magnetic resonance imaging data including T1-weighted, T2-weighted, and vascular imaging sequences, and raw image data from various modalities such as computed tomography angiography. During the reading process, the metadata header of the image file needs to be parsed to extract spatial positioning information such as pixel spacing, slice thickness, and image orientation. The images are then subjected to quality screening, for example, using a deep learning-based motion artifact detection model with a confidence threshold set to 0.95 to eliminate motion artifact interference and ensure that the spatial resolution of the images meets the accuracy requirements of three-dimensional sub-millimeter reconstruction with an isotropic voxel side length of no more than 0.5 mm, thereby establishing the individual anatomical geometric benchmark of the patient. Meanwhile, in acquiring dynamic physiological function data, continuous time-series signals are collected by connecting to data stream interfaces such as HL7 or IEEE11073PHD standard interfaces of clinical monitoring equipment or reading stored historical monitoring records. These signals include electrocardiogram signals with a sampling frequency of not less than 250Hz, electroencephalogram signals, invasive or non-invasive blood pressure waveforms with a sampling frequency of not less than 100Hz, hemodynamic indicators such as cardiac output and peripheral vascular resistance, and respiratory rate. The collected multi-channel signals are then processed for time-domain synchronization. Data sources with different sampling frequencies are aligned using a unified timestamp to eliminate clock drift between devices. The synchronization accuracy requirement is no more than 10 milliseconds to ensure the instantaneous consistency of boundary conditions during subsequent modeling.

[0075] Furthermore, when acquiring time-series metabolic parameter data, the patient's blood biochemical test results and body fluid metabolism analysis records are mainly extracted through the LIS system, including but not limited to blood glucose concentration, blood lactate level, blood oxygen saturation, acid-base balance, and electrolyte concentrations such as potassium and sodium ions. For discrete test data that are usually taken 1 to 4 times a day, a piecewise cubic Hermitian interpolation algorithm that ensures the continuity of the first derivative of the interpolation curve is used to fit it into a continuous metabolic curve that changes over time, so as to reflect the dynamic biochemical evolution of the internal environment. In addition, regarding the acquisition of molecular and genetic characteristic data, based on targeted gene testing reports such as whole-genome sequencing, exome sequencing, or next-generation sequencing with a patient coverage of at least 30-fold, bioinformatics analysis workflows such as using the GATK workflow to identify variants and annotate dbSNP, ClinVar, and PharmGKB databases are employed to identify specific single nucleotide polymorphisms, gene copy number variations, and the expression abundance of key genes such as FPKM values ​​from RNA-seq. In particular, genetic loci data that have been shown to be related to pharmacokinetics of CYP450 family genes, biomechanical properties of FBN1 genes related to connective tissue mechanical properties, and disease susceptibility are screened. Finally, in constructing the multi-source heterogeneous biomedical feature pool, standardized preprocessing operations were performed on the four types of heterogeneous data acquired above. These included removing environmental noise and outliers from the original signals, such as using a Butterworth bandpass filter of 0.5-40Hz for electrocardiogram signals, performing Z-score standardization on values ​​of different dimensions, and performing unified registration of spatiotemporal dimensions based on anatomical coordinates and absolute time axis. All dynamic data time axes were uniformly aligned to the zero point of the first image acquisition. The processed data were then structured and indexed according to anatomical labels such as "left ventricle" and "middle cerebral artery," physiological classifications such as "electrophysiology" and "hemodynamics," and time attributes. Finally, the data was aggregated into a dataset containing complete patient biomedical features that can be directly called by subsequent digital twin modeling algorithms and can be stored in HDF5 or Apache Parquet format.

[0076] Step S20: The static anatomical structure data is parsed to generate a three-dimensional anatomical structure model, and cross-scale fusion modeling is performed by combining the dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic characteristic data to form the initial parameter set of the personalized digital twin.

[0077] In this embodiment, the static anatomical structure data is parsed to generate a three-dimensional anatomical structure model, and the method is as follows:

[0078] Step S21: Spatial registration and voxel-level semantic fusion are performed on the multimodal static anatomical imaging data to generate unified three-dimensional volume data with tissue category annotations;

[0079] Step S22: Perform hierarchical segmentation on the three-dimensional volume data based on prior anatomical atlas, extract and reconstruct the boundary surfaces of different anatomical structures and substructures to form a three-dimensional mesh model that expresses the geometric morphology.

[0080] Step S23: Based on the three-dimensional mesh model, analyze and construct a multi-scale topological map describing the spatial adjacency relationship between organs and the physiological connectivity within tissues;

[0081] Step S24: The material property parameters characterizing the biophysical properties of the tissue are assigned to the corresponding units of the three-dimensional mesh model according to the tissue category label and the multi-scale topology map, thereby generating a computable three-dimensional anatomical structure model integrating geometric, topological and physical properties. The three-dimensional anatomical structure model provides the digital twin kernel engine with a discretized spatial domain for multi-physics coupling calculation.

[0082] In this embodiment, cross-scale fusion modeling is performed by combining the dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic characteristic data to form the initial parameter set of the personalized digital twin. The method is as follows:

[0083] Step S25: Construct a system state matrix based on the dynamic physiological function data and the time-series metabolic parameter data, and apply a system identification algorithm to solve the set of ordinary differential equations describing macroscopic physiological behavior to obtain a primary kinetic parameter set;

[0084] Step S26: Map the molecular and genetic feature data to a pre-set biological pathway knowledge graph. By calculating the centrality measure of specific biomolecular entities in the graph topology and the weights of their associated pathways, a set of genetic regulatory coefficients are generated.

[0085] Step S27: Using the genetic regulation coefficient, the primary dynamic parameter set is corrected through a preset parameter modulation function. The modulation function takes the genetic regulation coefficient as input and performs element-wise transformation on the state transition matrix of the primary dynamic parameter set to obtain the corrected system parameters under genetic background constraints.

[0086] Step S28: Establish a cross-scale mapping model with physical field smoothness and conservation law as joint optimization objectives, and map the corrected system parameters as the source field to the physical attribute target field of each discrete unit of the three-dimensional anatomical structure model to generate the initial parameter set of the personalized digital twin.

[0087] Specifically, in step S20, the static anatomical structure data is parsed to generate a three-dimensional anatomical structure model, and the dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic characteristic data are combined to perform cross-scale fusion modeling to form the initial parameter set of the personalized digital twin. The first step is to construct a high-fidelity geometric and physical base.

[0088] In step S21, the acquired multimodal static anatomical imaging data, specifically including computed tomography and magnetic resonance imaging data, are imported into the BRAINSFit module or SimpleElastix library of the 3DSlicer software, or other multimodal registration algorithm frameworks. This process first selects the modality with the highest resolution as the reference image, calculates the mutual information or normalized cross-correlation coefficient between it and other modal images, and iteratively optimizes the rigid or affine transformation matrix. The optimizer uses regular gradient descent with a maximum iteration count of 200 to achieve pixel-level strict alignment of different modal images in anatomical space, with a target registration error requirement of less than 1.0 mm. Subsequently, voxel-level semantic fusion is performed. For voxels at the same spatial location, a weighted fusion strategy is constructed based on the imaging advantages of different imaging modalities for soft and hard tissues. Specifically, a sigmoid weighting function based on grayscale thresholds is set. In high-density bone regions with CT values ​​greater than 200 Henle units, computed tomography (CT) values ​​are adopted with a weight of 0.8 to 1.0, while in soft tissue regions with CT values ​​between 0 and 100 Henle units, magnetic resonance imaging (MRI) values ​​are adopted with a high weight of, for example, 0.9. This generates a unified three-dimensional volume data with artifact removal, enhanced contrast, and each voxel carrying preliminary tissue probability information, with a uniform voxel size of 0.5 mm × 0.5 mm × 0.5 mm.

[0089] Next, step S22 is executed, employing a hierarchical segmentation method guided by prior anatomical atlases to process the aforementioned three-dimensional volume data. This method utilizes a B-spline freeform deformation algorithm with control point grid spacing set to 4 times the voxel size to map statistical shape models or probability maps containing standard anatomical structure annotations, obtained from public datasets such as ADNI or UKBIobank, or using internally constructed models, to the patient's individual space, providing initial spatial constraints for segmentation. Based on this, a fully convolutional neural network based on a U-Net architecture, such as using the nnU-Net framework for training and inference, or a level set evolution algorithm incorporating edge stopping functions, iterating 500 times with a time step of 0.1, is applied to refine the boundaries of organs, blood vessels, and lesions. After segmentation, a moving cube algorithm with isosurface thresholds set according to tissue type (e.g., setting the cardiac muscle to 50% maximum strength, or a Poisson surface reconstruction algorithm with a depth set to 10) is used to extract isosurfaces of various anatomical tissues and micro-substructures. The generated mesh is then iterated 5 times with a smoothing factor λ. smoothThe mesh extraction optimization, which uses a Laplacian smoothing of 0.5 and a quadratic error metric to reduce the number of facets by 50%, effectively removes stair-step artifacts and repairs non-manifold structures. The final result is a topologically correct 3D mesh model that can accurately represent the complex geometry of the patient's body, with tetrahedral mesh elements and an average Jacobian matrix determinant greater than 0.3.

[0090] In step S23, the inherent relationships within the 3D mesh model are further analyzed. A centerline extraction algorithm based on the fast traversal method is used to calculate the vessel centerline, extract the lumen cross-section, and calculate the distances between mesh patches. A contact threshold of 0.1 mm is set to identify the contact surfaces between organs, constructing an abstract multi-scale topological graph. This graph uses nodes to represent different organ cavities or functional units, and edges to represent blood circulation pathways, nerve conduction paths, or physical adjacency relationships. This transforms the complex geometric model into a mathematical atlas that describes the spatial adjacency relationships between organs and the physiological connectivity within tissues, and can be stored as a NetworkX graph object or a custom adjacency list, providing a logical connection foundation for subsequent fluid and biochemical simulations.

[0091] After completing the geometry and topology construction, step S24 is executed to give the model physical meaning, i.e., to generate a computable three-dimensional anatomical structure model. This process involves precisely mapping material property parameters characterizing the biophysical properties of tissues, including indicators such as Young's modulus, Poisson's ratio, thermal conductivity, and electrical conductivity, with benchmark values ​​obtained from biomechanical literature databases such as "Biomechanics of Living Tissues" or internal experimental data, onto each tetrahedral or hexahedral element of the three-dimensional mesh model, based on the tissue category labels generated in step S21 and the multi-scale topology map constructed in step S23. This mapping process is not a simple assignment, but rather establishes an empirical formula between grayscale values ​​and physical properties. For example, for bone elements, the formula is based on the Heinz units of CT values. The power-law relationship of Poisson's ratio ν=0.3 is used to calculate the non-uniform elastic modulus; E represents the elastic modulus of bone tissue. HU represents the Henle unit in computed tomography images.

[0092] For the vessel wall unit, anisotropic Neo-Hookean or Ogden hyperelastic material model parameters are assigned based on its diameter and spatial location. For example, the aortic wall parameters can be set as follows: Fiber orientation angular distribution parameters . The parameters represent the shear modulus of the vascular wall matrix material, k1 represents the stress parameters of collagen fibers, and k2 represents the dimensionless index parameters of collagen fibers.

[0093] The generated model is a digital entity that integrates geometric shape, topological logic, and nonlinear physical properties. It provides the digital twin kernel engine with a discretized spatial domain for multiphysics coupling calculations, which is ultimately exported in a simulation software compatible format such as FEBio or AbaqusINP.

[0094] After constructing the three-dimensional anatomical model and assigning its physical properties, the cross-scale fusion modeling steps S25 to S28 are further carried out. First, in step S25, dynamic physiological function data including heart rate and blood pressure, and time-series metabolic parameter data covering blood glucose and lactate levels are extracted to construct a state vector sequence describing the current operating state of the system, with a uniform sampling interval of 1 second. System identification techniques are applied, such as least squares, extended Kalman filtering, or neural network differential equation solvers. Specifically, a recursive least squares method with a forgetting factor is used. Set to 0.99 for a pre-defined zero-dimensional or one-dimensional lumped-parameter physiological model based on the Windkessel effect, for example, a differential equation of... The three-element RCR Windkessel model, where P is pressure, For inflow flow, For peripheral resistance, To ensure compliance, backfitting is performed to solve a set of ordinary differential equations that can reproduce macroscopic physiological behavior, along with a corresponding set of primary kinetic parameters for the systemic circulation model, which may include 15 to 30 state variables and corresponding parameters. These parameters represent the baseline physiological characteristics of the patient without considering microscopic genetic differences. Next, step S26 is executed to map the patient's molecular and genetic characteristic data to a pre-built biological pathway knowledge graph, sourced from KEGG, Reactome, or a custom medical knowledge base. This graph covers key biological processes such as signal transduction and metabolic circulation. The algorithm calculates metrics such as betweenness centrality and eigenvector centrality of specific biomolecules such as mutated genes and overexpressed proteins in the graph's topological network, and combines these with biological weights such as pathogenicity scores (ranging from 0 to 1) derived from literature for their respective pathways. This quantitatively assesses the potential impact of microscopic molecular changes on macroscopic function, generating a set of dimension-normalized genetic regulatory coefficients, each ranging from 0.5 to 2.0, with 1.0 representing no impact.

[0095] Subsequently, in step S27, this set of genetic regulatory coefficients is used to modify the primary kinetic parameter set. This is achieved through a preset parameter modulation function; for example, for pharmacokinetic parameters such as clearance rate CL, the modulation function is... ,in This is the genetic regulation coefficient, and α is an adjustable effect strength factor, which is set to 0.5 by default. This function defines an explicit mapping logic from microscopic genotype to macroscopic phenotypic parameters, using genetic regulatory coefficients as input variables to perform element-wise nonlinear transformations or linear weighting on the state transition matrix or reaction rate constant in the primary kinetic parameters. Specifically, the system employs a scaling factor adjustment mechanism. For example, if a loss-of-function mutation is detected in a drug-metabolizing enzyme gene such as CYP2C19, a decay coefficient less than 1, such as 0.3, is calculated based on its centrality score and directly multiplied into the kinetic parameters describing the drug clearance rate. This results in corrected system parameters that conform to real-time monitoring data while also incorporating genetic background constraints.

[0096] Finally, step S28 is executed to establish a cross-scale mapping model with the smoothness of the physical field and conservation laws, including mass and momentum conservation, as joint optimization objectives. This model acts as a bridge connecting zero-dimensional or one-dimensional lumped parameters with three-dimensional distributed parameters. Using the corrected system parameters as source field data, it maps them to the physical property target field of each discrete unit of the three-dimensional anatomical structure model through interpolation and Gaussian diffusion algorithms based on radial basis functions with Gaussian kernels and bandwidth parameters adaptively adjusted according to anatomical distance. This process ensures that macroscopic whole-body physiological parameters are rationally allocated to specific three-dimensional tissue spaces. For example, based on the principle of mass conservation, the whole-body cardiac output parameters are transformed into parabolic velocity distribution boundary conditions at specific blood vessel inlets, with peak velocity V... max It can be calculated using the formula (2 × cardiac output) / (π × vessel radius). 2 The calculation ultimately generates a structured dictionary or JSON object, where the key is the grid cell ID and the value is the personalized digital twin initial parameter set of the corresponding multi-physical attribute vector, completing a comprehensive integration from geometry to physics and from macro to micro.

[0097] Step S30: Drive the multi-physics field coupled digital twin kernel engine with the personalized digital twin initial parameter set to simulate the effect of preset intervention measures on the patient's full-parameter digital twin and output a set of predicted physiological state parameters.

[0098] In this embodiment, step S30 specifically includes the following steps:

[0099] Based on the mesh topology of the three-dimensional anatomical structure model, the blood rheological parameters and tissue mechanical modulus in the initial parameter set of the personalized digital twin are mapped to fluid domain units and solid domain nodes, respectively, to construct a multiphysics calculation model containing non-uniform material properties.

[0100] Step S31: Based on the multiphysics calculation model, establish a simultaneous solution mechanism including the Navier-Stokes equations, nonlinear constitutive equations, and reaction-convection-diffusion equations, and quantify the preset intervention measures as boundary conditions in the solution mechanism. Define the inlet concentration function or mass source term for drug treatment, and define the contact boundary or displacement constraint for physical intervention.

[0101] Step S311: Perform physical property-based region segmentation on the three-dimensional anatomical structure model to identify and mark the vascular lumen region, vascular wall region, and surrounding target tissue region.

[0102] Step S312: Based on the results of the region segmentation, assign region labels to the tetrahedral or hexahedral mesh units in the three-dimensional anatomical structure model, and assign specific material property values ​​to each mesh unit or node according to the parameter subset associated with each region label in the initial parameter set of the personalized digital twin. The property values ​​of the vascular lumen unit are set according to the blood rheological parameters, and the property values ​​of the vascular wall and target tissue units are set according to the tissue mechanical modulus.

[0103] Step S313: For mesh cells or nodes located at the junction of different material regions, use distance-based linear or nonlinear interpolation algorithms to smooth the material property values ​​and generate the multiphysics calculation model containing non-uniform material properties.

[0104] Step S32: Drive the simultaneous solution mechanism to perform fluid-structure interaction iteration within the time-separated walk, transferring the shear force and pressure on the fluid side to the solid side to calculate the deformation, and feeding back the deformation displacement to the fluid side to update the mesh, until the full-time-space numerical simulation is completed.

[0105] Step S321: A partitioned, strongly coupled iterative strategy is adopted to alternately solve the fluid dynamics and solid mechanics control equations and exchange boundary data in each global time step.

[0106] Step S322: Fix the current solid domain mesh shape, solve the Navier-Stokes equations to obtain the pressure field and velocity field of the fluid domain, and calculate the fluid load on the fluid-solid interface.

[0107] Step S323: Apply the fluid load to the corresponding boundary node of the solid domain and solve the nonlinear constitutive equation to obtain the deformation displacement field of the solid domain;

[0108] Step S324: Update the node coordinates of the fluid domain mesh using the arbitrary Lagrange-Euler method based on the deformation displacement field;

[0109] Step S325: Determine whether the forces and displacements on the fluid-structure interface meet the preset convergence conditions. If not, repeat the alternating solution based on the updated mesh and load until convergence, and then proceed to the next time step.

[0110] Step S33: For the physical field evolution data obtained from the simulation, a virtual sampling probe is set in the anatomical coordinate system, and the time series and statistical features of key variables are extracted through shape function interpolation, and the predicted physiological state parameter set is compiled and generated.

[0111] After obtaining a personalized initial parameter set for the digital twin containing detailed geometric and physical definitions, the system executes step S30. The core of this step is to use this parameter set to drive a multi-physics coupled digital twin kernel engine. This engine can be based on open-source libraries such as FEniCS and OpenFOAM, or a secondary development platform of commercial software ANSYS Fluent / Mechanical, to construct a highly realistic digital experimental environment. This process first maps the blood rheological parameters and tissue mechanical moduli from the parameter set to fluid domain elements and solid domain nodes, respectively, based on the mesh topology of the three-dimensional anatomical model, constructing a multi-physics computational model containing non-uniform material properties. Specifically, the system first executes step S31 and its sub-steps. In step S311, the three-dimensional anatomical model is divided into regions based on physical properties. By calculating the normal vector features of the mesh elements and combining them with the connected component labeling algorithm of the region growing algorithm, seed points are automatically selected based on tissue labels, and the similarity threshold is set to a material property difference of less than 10%. This automatically identifies and strictly distinguishes the vascular lumen region as the fluid computation domain, the vascular wall region as the solid deformation computation domain, and the surrounding target tissue region providing elastic support boundaries.

[0112] In step S312, the system assigns specific constitutive model parameters to each grid cell based on the region segmentation results. For fluid cells within the vascular lumen, to accurately describe the shear thinning characteristics of blood in microvessels or high-shear regions, the system abandons the simple Newtonian fluid assumption and instead adopts a Carreau-Yasuda or Casson non-Newtonian fluid constitutive model, such as the Carreau model. ,in, For apparent viscosity, Shear rate, The limiting viscosity at infinitely high shear rates. The limiting viscosity at zero shear rate. It is a time constant. It is the power-law exponent; its parameters can be set as follows: =0.0035 Pa·s, =0.056 Pa·s, =3.313s, n=0.3568. By setting specific parameters such as the limiting high shear viscosity, zero shear viscosity, and power-law exponent, a mathematical relationship between viscosity and shear rate dynamics is established. For the solid elements of the blood vessel wall and target tissue, the system adopts the Holzapfel-Gasser-Ogden (HGO) hyperelastic constitutive model, which describes the anisotropic mechanical behavior of soft tissue, and its strain energy function is: , For strain energy density, For matrix material parameters, For the first isocompressive strain invariant, These are the stress parameters for collagen fibers. This is a dimensionless index parameter for collagen fibers. The dispersion parameter represents the fiber orientation distribution. Let be the strain pseudo-invariant along the anisotropic direction.

[0113] Typical parameters include: C10 = 0.05 MPa, k1 = 0.2 MPa, k2 = 10, κ = 0.2. The fiber orientation is defined according to the local coordinate system. This model not only includes Mooney-Rivlin parameters describing matrix behavior but also defines a structural tensor characterizing the distribution direction of collagen fibers, thus accurately reflecting the nonlinear hardening behavior of the blood vessel wall under pulsating pressure. To eliminate the numerical singularity that may be caused by abrupt changes in stiffness at the interface between different material regions, the system executes step S313 to smooth the mesh node properties near the fluid-structure interface. An interpolation algorithm based on inverse distance weights is used to smooth three layers of elements on both sides of the interface. The weights are based on the formula... Calculate, where d is the distance to the interface. The constant is set to 0.001 mm to prevent division by zero, and a buffer layer for continuous physical property transition is established between different material property regions to ensure the numerical stability of the stiffness matrix assembly. Furthermore, in setting the boundary conditions, the system couples the three-element RCRWindkessel lumped parameter model, whose parameters are derived from the system identification results in step S25, through a user-defined function, such as near-end drag. Compliance C = 0.001 cm 5 / dyne, distal resistance Rd=1200dyne·s·cm -5 The resistance and compliance of the distal vascular bed are applied as implicit functions of pressure and flow to the outlet section of the three-dimensional model to reproduce the physiological pressure wave reflection effect.

[0114] After model construction and boundary definition are completed, the system executes step S32, driving the simultaneous solution mechanism to perform fluid-structure interaction iterative calculations within discrete time steps. Specifically, in this implementation, step S321 employs a partitioned strongly coupled iterative strategy based on arbitrary Lagrange-Euler descriptions, discretizing each cardiac cycle into 200 global time steps, i.e., time step size... The time step is approximately 5 milliseconds to balance computational efficiency and conservation accuracy. Within each time step, the system executes step S322, first freezing the current mesh displacement, discretizing the Navier-Stokes equations using the finite volume method, and decoupling the pressure and velocity fields using the PISO algorithm with pressure-velocity coupling and three internal iterations to solve the momentum and continuity equations. Then, it integrates to calculate the stresses, including pressure and viscous shear stresses, experienced by each surface element at the fluid-structure interface, according to the formula... Calculate the total traction force vector. Wherein, This represents the total traction force vector acting per unit area at the fluid-solid interface. For hydrostatic pressure, This represents the unit normal vector from the solid to the fluid at a point on the fluid-solid interface. The dynamic viscosity of the fluid. For the velocity gradient tensor, It is the transpose of the velocity gradient tensor.

[0115] Subsequently, the system executes step S323, transferring the total traction force vector as a Neumann boundary condition to the solid domain, and employing the finite element method combined with parameter β. N It is 0.25 and γ N The nonlinear dynamic equilibrium equations of the solid are solved using the Newmark-β time integration scheme with a time factor of 0.5. Where β... N and γ N These are the parameters of the integration algorithm, used to control the accuracy and stability of the numerical calculation. Given the large deformation characteristics of the blood vessel wall, the solution process includes Newton-Raphson iterations involving geometric and material nonlinearities, with the iteration tolerance residual set to 10 within each step. -3 The final output is the displacement field of the solid domain nodes.

[0116] To ensure the fluid mesh can adapt to the deformation of the solid without distortion, the system executes step S324, updating the fluid domain mesh based on the deformation displacement field of the solid. This step employs a mesh smoothing algorithm based on the diffusion equation to solve the equation. ,in, The mesh stiffness is set to 1 at the boundary and increased to 1000 internally to prevent excessive distortion. The boundary displacement is used as a source term to diffuse into the fluid domain. The Laplace equation is solved to smoothly redistribute the internal nodes while maintaining the mesh topology connectivity. For smooth displacement field The gradient;

[0117] Regarding the coupling convergence control at each time step, the system executes step S325 to determine whether the forces and displacements at the fluid-structure interface meet the preset convergence conditions. To accelerate convergence and prevent oscillations, the system introduces the Aitken dynamic relaxation algorithm, with a relaxation factor... = Dynamic calculation, where r is the residual. The dynamic relaxation factor is given in the k-th fluid-structure interaction iteration. The dynamic relaxation factor at the (k-1)th iteration. Let be the residual vector after the k-th iteration. This is the residual vector after the (k-1)th iteration.

[0118] By analyzing the relationship between the residual vectors of the previous and current iterations, the optimal relaxation factor is dynamically calculated to perform weighted updates on the interface displacements. The system continuously monitors the displacement and force residuals at the fluid-structure interface, and their relative residual norms are as follows: and until it falls below, for example, 10. -4 Only when a preset convergence threshold is set is the current time step considered converged and the process proceeds to the next time step, thus ensuring that the full-time-space numerical simulation strictly satisfies the laws of conservation of mass and momentum.

[0119] in, This represents the change in the predicted displacement of the fluid-structure interaction interface nodes between two consecutive fluid-structure interaction iterations. This indicates the displacement state of the interface at the start of the current global time step (i.e., before the start of the fluid-structure interaction iteration). This represents the difference in force vectors calculated at the interface between the fluid and the solid. This represents the interface force calculated at the start of the current global time step or after the first coupling iteration.

[0120] After the simulation completes the full-spatiotemporal evolution covering the preset cardiac cycle, for example, simulating three complete cardiac cycles and taking the data from the last cycle for steady-state analysis, step S33 is executed for post-processing. In this embodiment, the physical field evolution data obtained from the simulation is post-processed to generate the predicted physiological state parameter set. The method is as follows:

[0121] Within the anatomical coordinate system of the three-dimensional anatomical structure model, virtual sampling probes are set at preset key anatomical locations;

[0122] Using shape function interpolation techniques, the time series of key physical variables at the location of the virtual sampling probe are extracted from the discrete physical field evolution data output by the simultaneous solution mechanism;

[0123] Based on the extracted time series, advanced hemodynamic parameters including time-mean wall shear stress, oscillatory shear index, and relative residence time are calculated, wherein the relative residence time is calculated based on the functional relationship between the time-mean wall shear stress and the oscillatory shear index.

[0124] The time series of the key physical variables, the calculated advanced hemodynamic parameters, and the corresponding anatomical location coordinates and timestamps are compiled and integrated to generate a structured set of predicted physiological state parameters.

[0125] Specifically, at key locations within the anatomical coordinate system, such as the narrowest section at the stenosis, one diameter downstream of the bifurcation point, the aneurysm neck, and the aneurysm top, virtual sampling probes defined as "PointProbe" in the simulation software or defined by coordinates are set. Shape function interpolation techniques are used to extract time series data from discrete field data that match the output frequency, such as 200 points per second. Further advanced hemodynamic parameters are calculated, including time-averaged wall shear stress in Pascals, oscillatory shear index ranging from 0 to 0.5, and parameters calculated according to formulas. The calculated relative dwell time. Wherein, Indicates relative stay time. This represents the average magnitude of the wall shear stress vector over a complete cardiac cycle. It is a dimensionless index representing the degree of oscillation of the direction of wall shear stress within a cardiac cycle.

[0126] These derived parameters, along with basic pressure, flow velocity, and deformation data, are compiled into a set of predicted physiological state parameters that can be stored as CSV files or database records containing timestamps, location coordinates, and various parameters. This provides a quantitative assessment basis that includes flow field details and biomechanical responses for subsequent intervention program optimization.

[0127] Step S40: Construct an inverse optimization problem based on the difference between the preset clinical goal and the predicted physiological state parameter set, and introduce dynamic physiological state constraints based on real-time monitoring data of the patient during the optimization process to solve and generate optimized intervention plan parameters.

[0128] In this embodiment, the method of step S40 includes:

[0129] Step S41: Based on the quantitative difference between the predicted physiological state parameter set and the preset clinical target, construct an inverse optimization problem with minimizing the difference as the objective function;

[0130] Step S42: The real-time monitoring data of the patient is converted into a feasible domain boundary condition that evolves over time, and this boundary condition is introduced as an inequality constraint term into the inverse optimization problem to form a dynamic constraint optimization model.

[0131] Step S43: The dynamic constraint optimization model is iteratively solved using a gradient-based optimization algorithm. In each iteration, the feasible region boundary conditions are updated synchronously to reflect the latest monitoring data, and the search direction is guided by calculating the constraint violation penalty term until the algorithm converges.

[0132] Step S44: Output the decision variable values ​​corresponding to the convergence, as parameters of the optimized intervention scheme.

[0133] Step S42 includes:

[0134] Step S421: The real-time monitoring data of the patient is processed in real time, and the trend characteristics and fluctuation range of key physiological parameters are extracted by sliding time window;

[0135] Step S422: Based on the extracted trend features and fluctuation range, establish a dynamic feasible region describing the safety boundary of the physiological state, and model the dynamic feasible region as a time-varying inequality constraint of the decision variables in the inverse optimization problem;

[0136] Step S423: During the iterative solution process of the optimization algorithm, the specific parameters of the time-varying inequality constraint are updated according to the latest window of the monitoring data to ensure that the parameters of the generated optimized intervention plan always meet the patient's current physiological state boundary.

[0137] After obtaining a set of predicted physiological state parameters containing high-fidelity flow field details and biomechanical responses, the system proceeds to step S40. The core of this step lies in transforming the predictive capabilities of the digital twin into clinical decision support capabilities. The system does not simply select the optimal solution, but rather constructs a complex inverse optimization problem incorporating a real-time feedback mechanism. Specifically, in step S41, the system first defines a quantified objective function for optimal selection. This function is constructed based on the difference between the predicted set of physiological state parameters and a preset clinical target. The preset clinical target could specifically be restoring the wall shear stress of the narrowed vessel segment to a normal physiological range, such as 1.5 ± 0.5 Pa, or maintaining the drug concentration in the target organ within a therapeutic window, such as 2.0 ± 0.3 μg / mL. To comprehensively consider multi-dimensional treatment effects, this objective function is typically designed in a weighted least squares form, i.e. ,in Decision variables include, for example, stent diameter and drug dosage. It is the i-th prediction parameter. It is the i-th preset clinical target value. The weights are determined based on the priority of clinical indications. For example, in aneurysm treatment, reducing the oscillatory shear index with a weight w1 of 0.6 will have a higher weight than improving the mean intraneural flow velocity with a weight w2 of 0.4. The system aims to find a set of intervention parameters that best approximates the output of the digital twin model to the ideal clinical state by minimizing this objective function.

[0138] Unlike traditional static optimization, this implementation introduces a dynamic constraint mechanism based on real-time patient monitoring data in step S42, which is particularly crucial for intensive care or intraoperative navigation scenarios. This mechanism first executes step S421, where the system connects to the patient's bedside monitoring equipment or wearable sensor data stream. It then uses a sliding time window with a specific width of 30 seconds to 1 minute to sample and process continuous physiological signals such as heart rate, blood pressure, and blood oxygen saturation in real time. Within the window, the system applies adaptive filtering algorithms such as Kalman filtering or wavelet denoising to remove signal noise and extracts the immediate trend characteristics of key physiological parameters, such as the average rate of change (slope) over 5 seconds, and the fluctuation range represented by the standard deviation, thereby obtaining a dynamic profile reflecting the patient's current physiological stability and reserve capacity.

[0139] Then, in step S422, the system maps the extracted trend and fluctuation features to geometric constraint boundaries in the decision variable space, establishing a dynamic feasible region describing the safety boundary of the physiological state. For example, if monitoring data shows that the patient's systolic blood pressure is decreasing and the fluctuation is increasing, such as a decrease in mean blood pressure of more than 10 mmHg and an increase in standard deviation of more than 5 mmHg within 10 minutes, the system will automatically reduce the allowable upper limit of the dosage of vasodilators, lowering the original upper limit D. max The parameters can be dynamically adjusted from 10 mg / h to 7 mg / h, or the range of certain high-risk values ​​for interventional surgical parameters can be limited, such as reducing the upper limit of stent expansion pressure from 16 atmospheres to 14 atmospheres. Mathematically, this is modeled as a time-varying inequality constraint in an inverse optimization problem, i.e., g... j(u,t) ≤0, j=1,2,...m, where t represents the actual time during optimization, ensuring that the search space is always limited within the patient's current physiological tolerance limit, g j(u,t) Let be the j-th time-varying inequality constraint function defined on the decision variable space u at time t, where m is the total number of constraints.

[0140] After constructing the objective function and dynamic constraints, the system executes step S43, using a gradient-based optimization algorithm to iteratively solve the dynamic constraint optimization model. The method is as follows:

[0141] In each optimization iteration, the gradient of the objective function with respect to the decision variables is first calculated to determine the search direction;

[0142] At the same time, based on the latest real-time monitoring data window, the specific parameters of the time-varying inequality constraints are updated to dynamically adjust and optimize the feasible region;

[0143] During the solution process, a penalty function mechanism or augmented Lagrange multiplier method is used to handle constraint violations. The constraint violation amount is added as a penalty term to the objective function, or the search path is guided by multiplier updates to avoid constraint violation areas caused by changes in the patient's physiological state.

[0144] The above process is executed iteratively until the gradient norm of the objective function is less than the preset convergence tolerance, and all dynamic constraints meet the optimization convergence conditions. Then the algorithm is considered to have converged.

[0145] Specifically, the nonlinear dynamically constrained optimization model is solved iteratively using efficient gradient-based optimization algorithms, such as sequential quadratic programming implemented through NLopt or SciPy libraries, or adjoint equation-based methods that automatically calculate gradients using libraries like dolfin-adjoint. In each optimization iteration, the system not only calculates the gradient of the objective function with respect to the decision variables using finite difference methods or automatic differentiation to determine the search direction, but more importantly, it executes step S423 to simultaneously check the latest real-time monitoring data window and update the specific parameters of the time-varying inequality constraints based on the latest physiological state. For example, after each digital twin simulation, which takes approximately several minutes, the latest 30 seconds of monitoring data is read from the data buffer and the constraint boundaries are recalculated. This means that the feasible region of the optimization algorithm is not fixed but dynamically adjusted in real time as the patient's state evolves. To handle this dynamism, the system integrates a penalty function mechanism into the optimization solver, adding constraint violations to the objective function. ,in It is an initial penalty factor of 1 that is multiplied by 10 after each violation, or the augmented Lagrange multiplier method, which strongly guides the search path to avoid the taboo areas newly generated due to the deterioration of the patient's condition by calculating the constraint violation penalty term. For the decision variable vector, The original objective function is... Let j be the constraint function for the j-th inequality.

[0146] This process is repeated until the gradient norm of the objective function is less than, for example, 10. -5If the preset tolerance is met and all dynamic constraints are satisfied (i.e., the Karush-Kuhn-Tucker (KKT) conditions are met), the algorithm is considered converged. At this point, the system executes step S44 and outputs the decision variable values ​​corresponding to the convergence time. These values ​​are the parameters of the optimized intervention plan under the premise of balancing ideal treatment effect and real-time physiological safety of the patient. Specifically, they cover key variables such as the optimal drug infusion rate of 5.2 mg / h, and the optimal stent implantation size and position of 4.0 mm stent diameter with the proximal positioning marker 3.1 mm from the branch port. This completes the closed-loop operation from simulation prediction to intelligent planning.

[0147] Step S50: Input the optimized intervention plan parameters into the digital twin kernel engine for verification simulation, and generate a treatment simulation report based on the set of verification physiological state parameters output by the simulation.

[0148] In this embodiment, step S50 specifically includes:

[0149] Step S51: Configure the parameters of the optimized intervention scheme as the input boundary conditions or source terms of the digital twin kernel engine, drive the kernel engine to perform a verification numerical simulation based on the same physical mechanism as the original simulation but with optimized parameters, and output a set of verification physiological state parameters containing complete spatiotemporal evolution data.

[0150] Step S52 involves post-processing analysis of the verified physiological state parameter set, including calculating its achievement degree with the preset clinical goal, the improvement magnitude with the predicted physiological state parameter set, and assessing its long-term evolution trend and potential risks based on the kinetic stability criterion.

[0151] Step S521: Quantify the difference between the set of verified physiological state parameters and the preset clinical target. Specifically, calculate the difference norm between the set of verified physiological state parameters and the preset clinical target vector, and calculate the target achievement rate based on the ratio of the difference norm to the preset clinical target vector norm.

[0152] Step S522: Retrieve the predicted physiological state parameter set as a baseline, and for the same key physiological indicators, calculate the relative change of the values ​​in the verification physiological state parameter set relative to the baseline values. This relative change represents the improvement brought about by the optimized intervention plan.

[0153] Step S523: Based on the system evolution data contained in the set of verified physiological state parameters, the long-term stability of the physiological system after intervention is evaluated by constructing the system dynamic state trajectory or analyzing the eigenvalues ​​of the Jacobian matrix of the linearized system model, and potential unstable modes whose real parts of the eigenvalues ​​exceed the preset stability threshold are identified and marked.

[0154] Step S53: The results of the post-processing analysis, together with the parameters of the optimized intervention plan and the visualization data of the key simulation process, are structurally integrated according to a preset template to generate the treatment simulation report containing quantitative assessment, plan comparison and risk warning information.

[0155] To ensure that the final treatment plan is not only mathematically optimal but also completely rational and safe in terms of biophysical mechanisms, the system does not directly output the optimized intervention parameters to the doctor after obtaining them. Instead, it enters a crucial closed-loop verification phase, namely, step S50. This process first executes step S51, where the system remaps the mathematically converged optimized intervention parameters calculated in the previous step—such as the specific drug infusion time curve defined in CSV file format, the exact coordinates and expansion diameter of the stent implantation with modified geometric parameters and boundary conditions in the simulation software—back to physical space and configures them as input boundary conditions or mass source terms for the digital twin kernel engine. At this point, the system restarts the multiphysics coupling engine, driving the kernel engine to perform a high-precision verification numerical simulation.

[0156] It is worth noting that this simulation strictly follows the same physical mechanism as the initial state assessment, namely the same governing equations, the same material constitutive model, the same fluid-structure interaction strategy, and the same convergence criterion. The only variable is the input intervention parameters. This step aims to eliminate approximation errors that may be introduced by the optimization algorithm. By solving the full-order Navier-Stokes equations and nonlinear solid constitutive equations, a set of verification physiological state parameters is generated, which includes the complete spatiotemporal evolution of all physical quantities such as flow velocity, pressure, shear stress, and tube wall deformation throughout the entire cardiac cycle. The data format is the same as the output of step S33.

[0157] After acquiring detailed validation data, the system proceeds to step S52, performing in-depth post-processing analysis on the massive four-dimensional spatiotemporal data (three-dimensional space plus one-dimensional time). The system first calculates the physiological parameters obtained through validation simulations against preset clinical targets, such as the target wall shear stress value (TAWSS). target =1.5 Pa, target fractional flow reserve (FFR) target A similarity score greater than 0.80 is used to quantify the Euclidean distance or correlation coefficient between the two. The degree of achievement can be determined using the formula... The calculations are performed, and the results are presented as a percentage of the target achievement rate. To verify the physiological state parameter vector, This is a pre-defined clinical target vector.

[0158] Simultaneously, the system retrieves the original prediction data before the intervention was applied, calculates the relative improvement of key hemodynamic parameters, and the improvement is based on the formula... Calculate, where, Baseline status parameter values.

[0159] For example, calculating the percentage reduction in the oscillatory shear index or the multiple increase in local drug concentration. More importantly, step S52 also introduces a kinetic stability criterion. By calculating the phase diagram trajectory of the flow field evolution, such as drawing a velocity-pressure phase diagram at key points, or evaluating the eigenvalue distribution of the system's Jacobian matrix, the system is linearized through small perturbation analysis, and the real parts of the eigenvalues ​​are calculated to determine whether there are positive real parts leading to instability. This allows for the assessment of whether there is a potential risk of flutter, increased turbulence, or aneurysm wall resonance in the vascular system after treatment under long-cycle physiological pulsation, which can be extended to 10 to 20 cardiac cycles. This provides a prediction of the long-term stability and safety of the treatment effect, and marks any unstable modes with real parts of eigenvalues ​​greater than -0.1.

[0160] Finally, the system executes step S53, transforming the complex quantitative analysis results into intuitive and easy-to-read clinical decision support documents. The system utilizes volume rendering and streamline tracing techniques, employing ParaView, VTK, or Matplotlib for visualization, to convert key simulation processes into visual images or animations. Examples include MP4 videos comparing time-averaged wall shear stress cloud maps before and after optimization, and GIFs depicting changes in intratumoral streamlines. These visually demonstrate changes in blood flow traces, the elimination process of high wall shear stress regions, and drug diffusion distribution cloud maps within the tissue. The system uses pre-defined standardized templates, such as Jinja2-based HTML or Word document templates, to integrate patient information, histograms comparing key indicators before and after optimization, risk warning signals (e.g., low safety margin areas highlighted in red), and accompanying text descriptions such as "The oscillatory shear index at the aneurysm neck remains above 0.3, indicating a long-term risk of thrombosis," along with the final recommended intervention parameters, such as "It is recommended to use a 3.5mm × 18mm drug-eluting stent, positioned 15 mm from the opening of the anterior descending artery, with a release pressure of 12 atmospheres. It is expected that the postoperative fractional flow reserve at the stenosis will increase to 0.88, and the time-averaged wall shear stress will recover to 1.6 Pa." This process automatically compiles and generates a comprehensive treatment simulation report containing quantitative assessment conclusions, cross-sectional comparisons of multiple treatment options, and potential risk warnings. This report can be output in PDF, HTML, or the DICOM-SR (structured report) format, directly embedded in the hospital's image archiving and communication system, providing clinicians with scientifically validated data through digital twins to ultimately determine surgical or medication plans.

[0161] The following is a specific example, which uses a virtual patient with coronary artery stenosis to fully implement the method described in this invention, aiming to plan the optimal percutaneous coronary intervention for the patient.

[0162] Those skilled in the art will understand that the following examples are merely illustrative of specific applications of the present invention and not intended to limit its scope of protection; the specific numerical values, software names, and model parameters involved are all exemplary.

[0163] Optimization simulation of PCI procedure for patients with coronary artery stenosis

[0164] A 65-year-old male patient presented with exertional angina and was initially diagnosed with suspected coronary atherosclerotic heart disease. To accurately assess his condition and plan a personalized treatment strategy, the method and system described in this invention were introduced.

[0165] Step S10 involves comprehensive data acquisition and integration. The system acquired high-resolution coronary computed tomography (CT) images of the patient with a slice thickness of 0.5 mm from the hospital's PACS / HIS system, serving as static anatomical data. Simultaneously, monitoring equipment was connected to collect Holter monitoring records and non-invasive continuous blood pressure monitoring data, serving as dynamic physiological function data. Recent blood test reports were reviewed to obtain time-series metabolic parameters such as blood lipids and blood glucose. Furthermore, targeted gene sequencing was performed on the patient, yielding molecular and genetic characteristic data, including key genes affecting the metabolism of antiplatelet drugs such as clopidogrel, such as CYP2C19. All data, after timestamp alignment and format standardization, were integrated into a multi-source heterogeneous biomedical feature pool indexed by the patient ID.

[0166] Step S20 is executed to construct a personalized digital twin model. Using CCTA data, a three-dimensional geometric model of the patient's left coronary artery tree is accurately reconstructed in medical image processing software such as MimicsInnovation Suite through threshold segmentation and manual correction. Specifically, a mixed plaque approximately 15 mm long, causing 75% luminal stenosis, located in the mid-segment of the left anterior descending artery (LAD), is clearly delineated. This geometric model is then imported into finite element preprocessing software such as AnsysSpaceClaim & Meshing, generating a high-quality computational mesh containing approximately 1.5 million tetrahedral elements, and marking the blood fluid domain, vessel wall solid domain, and plaque solid domain. Subsequently, dynamic physiological data (heart rate, blood pressure) are used to calibrate a three-element Windkessel model, which serves as the physiological boundary condition for the distal coronary artery exit. Its parameters are obtained through least-squares fitting, specifically including proximal resistance Rp equal to 98 dyne·s·cm. -5 The compliance is 0.0011 cm. 5 / dyne, and the distal resistance Rd equals 1150 dyne·s·cm -5These parameters collectively constitute the initial set of kinetic parameters. Crucially, systematic analysis revealed that the patient carries a loss-of-function CYP2C19 allele, leading to a significantly reduced ability to metabolize the antiplatelet drug clopidogrel. The system mapped this genetic trait to a pharmacokinetic (PK / PD) knowledge graph, generating a drug response regulation coefficient set to 0.4, implying only 40% of the standard drug response. This coefficient was used to revise the drug release and tissue uptake equations in the drug-eluting stent (DES) model, resulting in revised system parameters reflecting the genetic background. Finally, these parameters were mapped to a 3D model, assigning Carreau non-Newtonian fluid properties to the blood and an anisotropic HGO hyperelastic material model to the vessel wall, forming a personalized digital twin initial parameter set encompassing complete anatomical, physiological, and genetic characteristics.

[0167] Step S30 is executed to perform a pre-treatment simulation under baseline conditions. The initial parameter set is loaded into a multiphysics coupled simulation engine, such as the coupled solver of ANSYS Fluent and Mechanical, to simulate coronary hemodynamics during a complete cardiac cycle at a heart rate of 68 beats / min. The simulation results, i.e., the predicted physiological state parameter set, show that the calculated fractional flow reserve (FFR) distal to the stenosis is 0.72, lower than the clinical intervention threshold of 0.80, confirming the hemodynamic significance of the stenosis. Further analysis revealed that the wall shear stress (WSS) at the stenotic laryngeal region was as high as 9.2 Pa, while a large area of ​​low WSS (below 0.5 Pa) and a high-risk area with an oscillating shear index (OSI) above 0.3 appeared in the post-stenotic region. These are key risk factors leading to plaque progression and thrombosis.

[0168] Next, we proceed to the core step S40, which involves dynamic constraint optimization of the treatment plan. The preset clinical objectives are: postoperative FFR > 0.90 and minimization of the low WSS region area. The decision variables are set as three key parameters of the stent to be implanted during PCI: diameter D, ranging from 3.0 to 4.0 mm; length L, ranging from 18 to 28 mm; and longitudinal placement P. The system is designed to minimize... This is the inverse optimization problem of the objective function, where w are the weight coefficients. It is a key hemodynamic indicator, specifically referring to the total area of ​​all regions on the simulated vessel wall where the wall shear stress (WSS) is below a preset clinical threshold. During the optimization iteration process, the system incorporated real-time vital sign monitoring data from the patient under simulated surgical conditions. In one iteration, when the optimization algorithm attempted a more aggressive stent size, the patient's real-time monitoring data showed a compensatory increase in heart rate and a slight decrease in blood pressure. The system immediately converted this information into a dynamic constraint: in subsequent optimization searches, any stent parameter combination that causes a decrease in simulated mean arterial pressure exceeding 15% will be penalized or excluded. Using a sequential quadratic programming algorithm, after considering this dynamic physiological safety boundary, the algorithm converged after approximately 200 iterations, generating the optimal intervention parameters: a drug-eluting stent with a diameter of 3.5 mm and a length of 24 mm was selected, and its proximal marker should be precisely located 48.2 mm from the LAD ostium.

[0169] Step S50 is executed to perform a verification simulation and generate a report. The system applies the optimized stent parameters to the digital twin model, updates the geometric model by simulating the virtual stent implantation process, and then drives the kernel engine to perform a complete verification simulation. The set of verification physiological parameters output by the simulation shows that after implanting the stent with optimized parameters, the FFR value successfully increased to 0.91; the low WSS area in the post-stenotic region decreased by 95%, and the blood flow pattern was significantly improved. The system then automatically generates a graphic treatment simulation report. The report includes: a 3D visualization comparison of the coronary artery model before and after optimization, significant changes in the WSS cloud map, a quantitative improvement in the FFR value from 0.72 to 0.91, and clearly indicates the recommended stent type, size, and precise implantation coordinates. The report also includes a risk assessment: based on the optimized hemodynamic environment, the long-term risk of in-stent restenosis is predicted to be low. This report provided the cardiology team with strong, quantifiable decision support, enabling them to clearly anticipate the advantages, disadvantages, and potential risks of different treatment options before the actual surgery, thus allowing them to select the safest and most effective personalized treatment strategy for the patient.

[0170] See Figure 2 A system for constructing a full-parameter digital twin of a patient and performing treatment simulation according to a second embodiment of the present invention is used to execute a method for constructing a full-parameter digital twin of a patient and performing treatment simulation. The system includes:

[0171] The data integration module is configured to acquire and integrate static anatomical structure data, dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic characteristic data of the target patient to construct a multi-source heterogeneous biomedical feature pool.

[0172] The digital twin modeling module is configured to parse the static anatomical structure data to generate a three-dimensional anatomical structure model, and combine the dynamic physiological function data, time-series metabolic parameter data and molecular and genetic feature data to perform cross-scale fusion modeling to form the personalized digital twin initial parameter set.

[0173] The multi-physics coupling simulation module is configured to drive the multi-physics coupling digital twin kernel engine with the personalized digital twin initial parameter set, simulate the effect of preset intervention measures on the patient's full-parameter digital twin, and output a set of predicted physiological state parameters.

[0174] The dynamic constraint optimization module is configured to construct an inverse optimization problem based on the difference between the preset clinical goal and the predicted physiological state parameter set, and introduce dynamic physiological state constraints based on real-time monitoring data of patients during the optimization process to solve and generate optimized intervention plan parameters.

[0175] The verification and report generation module is configured to input the parameters of the optimized intervention plan into the digital twin kernel engine for verification simulation, and generate a treatment simulation report based on the set of verification physiological state parameters output by the simulation.

[0176] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working process and related explanations of the methods described above can be found in the corresponding processes in the foregoing system embodiments, and will not be repeated here.

[0177] It should be noted that the system for constructing a full-parameter digital twin of a patient and simulating treatment provided in the above embodiments is only an example of the division of the above functional modules. In practical applications, the above functions can be assigned to different functional modules as needed, that is, the modules or steps in the embodiments of the present invention can be further decomposed or combined. For example, the modules in the above embodiments can be merged into one module, or further divided into multiple sub-modules to complete all or part of the functions described above. The names of the modules and steps involved in the embodiments of the present invention are only for distinguishing the various modules or steps and are not considered as an improper limitation of the present invention.

[0178] A device according to a third embodiment of the present invention includes:

[0179] At least one processor;

[0180] and a memory communicatively connected to at least one of the processors;

[0181] The memory stores instructions that can be executed by the processor to implement the above-described method for constructing a full-parameter digital twin of a patient and performing treatment simulation.

[0182] A fourth embodiment of the present invention provides a computer-readable storage medium storing computer instructions, which are executed by the computer to implement the above-described method for constructing a full-parameter digital twin of a patient and performing treatment simulation.

[0183] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working process and related descriptions of the storage device and processing device described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0184] The terms “first”, “second”, etc., are used to distinguish similar objects, not to describe or indicate a specific order or sequence.

[0185] The term "comprising" or any other similar term is intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus / device that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent in such process, method, article, or apparatus / device.

[0186] The technical solution of the present invention has been described above with reference to the preferred embodiments shown in the accompanying drawings. However, it will be readily understood by those skilled in the art that the scope of protection of the present invention is obviously not limited to these specific embodiments. Without departing from the principles of the present invention, those skilled in the art can make equivalent changes or substitutions to the relevant technical features, and the technical solutions after these changes or substitutions will all fall within the scope of protection of the present invention.

Claims

1. A method of building a patient full parameter digital twin and performing therapy simulation, characterized in that, The method includes: Acquire and integrate static anatomical data, dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic characteristic data of target patients to construct a multi-source heterogeneous biomedical feature pool; The static anatomical structure data is analyzed to generate a three-dimensional anatomical structure model, and combined with the dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic feature data to perform cross-scale fusion modeling, forming a personalized digital twin initial parameter set; The personalized digital twin initial parameter set drives a multi-physics coupled digital twin kernel engine to simulate the effect of preset intervention measures on the patient's full-parameter digital twin, and outputs a set of predicted physiological state parameters: Based on the mesh topology of the three-dimensional anatomical structure model, the blood rheological parameters and tissue mechanical modulus in the initial parameter set of the personalized digital twin are mapped to fluid domain units and solid domain nodes, respectively, to construct a multiphysics calculation model containing non-uniform material properties. Based on the multiphysics calculation model, a simultaneous solution mechanism including the Navier-Stokes equations, nonlinear constitutive equations, and reaction-convection-diffusion equations is established. The preset intervention measures are quantified as boundary conditions in the solution mechanism. For drug treatment, an inlet concentration function or mass source term is defined, and for physical intervention, a contact boundary or displacement constraint is defined. The simultaneous solution mechanism is driven to perform fluid-structure interaction iteration within the time dispersive walk, transferring the shear force and pressure on the fluid side to the solid side to calculate the deformation, and feeding back the deformation displacement to the fluid side to update the mesh, until the full spatiotemporal numerical simulation is completed; For the physical field evolution data obtained from the simulation, a virtual sampling probe is set in the anatomical coordinate system, and the time series and statistical features of key variables are extracted by shape function interpolation, and the predicted physiological state parameter set is compiled and generated. An inverse optimization problem is constructed based on the difference between the preset clinical goals and the predicted physiological state parameter set. Dynamic physiological state constraints based on real-time patient monitoring data are introduced during the optimization process to solve and generate optimized intervention plan parameters. The optimized intervention parameters are input into the digital twin kernel engine for verification simulation, and a treatment simulation report is generated based on the set of verification physiological state parameters output by the simulation.

2. The method of building a patient's full parameter digital twin and conducting therapy simulation of claim 1, wherein, The method for parsing the static anatomical data to generate a three-dimensional anatomical model is as follows: Spatial registration and voxel-level semantic fusion of multimodal static anatomical imaging data are performed to generate unified 3D volume data with tissue category annotations. The three-dimensional volume data is subjected to hierarchical segmentation based on prior anatomical atlases, and the boundary surfaces of different anatomical structures and substructures are extracted and reconstructed to form a three-dimensional mesh model that expresses the geometric morphology. Based on the aforementioned three-dimensional mesh model, a multi-scale topological map describing the spatial adjacency relationships between organs and the physiological connectivity within tissues is analyzed and constructed. The material property parameters characterizing the biophysical properties of the tissue are assigned to the corresponding units of the three-dimensional mesh model according to the tissue category label and the multi-scale topology map, thereby generating a computable three-dimensional anatomical structure model that integrates geometric, topological and physical properties. The three-dimensional anatomical structure model provides the digital twin kernel engine with a discretized spatial domain for multi-physics coupling calculation.

3. The method of building a patient's full parameter digital twin and conducting therapy simulation of claim 2, wherein, By combining the dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic characteristic data to perform cross-scale fusion modeling, the initial parameter set for the personalized digital twin is formed. The method is as follows: Based on the dynamic physiological function data and the time-series metabolic parameter data, a system state matrix is ​​constructed, and a system identification algorithm is applied to solve the set of ordinary differential equations describing macroscopic physiological behavior to obtain a primary set of kinetic parameters. The molecular and genetic feature data are mapped to a pre-set biological pathway knowledge graph. By calculating the centrality measure of specific biomolecular entities in the graph topology and the weights of their associated pathways, a set of genetic regulatory coefficients are generated. Using the genetic regulation coefficient, the primary dynamic parameter set is corrected through a preset parameter modulation function. The parameter modulation function takes the genetic regulation coefficient as input and performs element-wise transformation on the state transition matrix of the primary dynamic parameter set to obtain the corrected system parameters under genetic background constraints. A cross-scale mapping model is established with physical field smoothness and conservation law as joint optimization objectives. The modified system parameters are used as the source field and mapped to the physical attribute target field of each discrete unit of the three-dimensional anatomical structure model to generate the initial parameter set of the personalized digital twin.

4. The method of building a patient's full parameter digital twin and conducting therapy simulation of claim 1, wherein, The multiphysics computation model is constructed as follows: The three-dimensional anatomical model is segmented based on physical properties to identify and mark the vascular lumen region, vascular wall region, and surrounding target tissue region. Based on the results of the region segmentation, tetrahedral or hexahedral mesh units in the three-dimensional anatomical structure model are assigned region labels, and based on the parameter subset associated with each region label in the personalized digital twin initial parameter set, specific material property values ​​are assigned to each mesh unit or node, wherein the property values ​​of the vascular lumen unit are set according to the blood rheological parameters, and the property values ​​of the vascular wall and target tissue units are set according to the tissue mechanical modulus. For mesh cells or nodes located at the boundaries of different material regions, distance-based linear or nonlinear interpolation algorithms are used to smooth material property values ​​and generate a multiphysics computation model that includes non-uniform material properties.

5. The method of building a patient's full parameter digital twin and conducting therapy simulation of claim 1, wherein, The method to drive the simultaneous solution mechanism to perform bidirectional fluid-structure interaction iteration within the time dispersive walk is as follows: A partitioned, strongly coupled iterative strategy is adopted to alternately solve the fluid dynamics and solid mechanics control equations and exchange boundary data in each global time step. With the current solid domain mesh shape fixed, solve the Navier-Stokes equations to obtain the pressure and velocity fields of the fluid domain, and calculate the fluid loads at the fluid-solid interface; The fluid load is applied to the corresponding boundary node of the solid domain, and the deformation displacement field of the solid domain is obtained by solving the nonlinear constitutive equation. Based on the deformation displacement field, the node coordinates of the fluid domain mesh are updated using an arbitrary Lagrange-Euler method; Determine whether the forces and displacements at the fluid-structure interface meet the preset convergence conditions. If not, perform the alternating solution again based on the updated mesh and load until convergence, and then proceed to the next time step.

6. The method of building a patient's full parameter digital twin and conducting therapy simulation of claim 1, wherein, An inverse optimization problem is constructed based on the difference between the preset clinical goals and the predicted physiological state parameter set. Dynamic physiological state constraints based on real-time patient monitoring data are introduced during the optimization process to solve for and generate optimized intervention plan parameters. The method is as follows: Based on the quantitative difference between the predicted physiological state parameter set and the preset clinical target, an inverse optimization problem is constructed with minimizing the difference as the objective function; The real-time patient monitoring data is converted into a feasible domain boundary condition that evolves over time, and the feasible domain boundary condition is introduced as an inequality constraint term into the inverse optimization problem to form a dynamic constraint optimization model. The dynamic constraint optimization model is solved iteratively using a gradient-based optimization algorithm. In each iteration, the boundary conditions of the feasible region are updated synchronously to reflect the latest monitoring data. The search direction is guided by calculating the constraint violation penalty term until the algorithm converges. The decision variable values ​​corresponding to the convergence are output as parameters of the optimized intervention scheme.

7. The method of building a patient's full parameter digital twin and conducting therapy simulation of claim 6, wherein, The real-time patient monitoring data is transformed into a feasible domain boundary condition that evolves over time, and this boundary condition is introduced as an inequality constraint term into the inverse optimization problem. The method is as follows: The patient's real-time monitoring data is processed in real time, and the trend characteristics and fluctuation range of key physiological parameters are extracted through a sliding time window; Based on the extracted trend features and fluctuation range, a dynamic feasible region describing the safety boundary of physiological state is established, and the dynamic feasible region is modeled as a time-varying inequality constraint of decision variables in the inverse optimization problem. During the iterative solution process of the optimization algorithm, the specific parameters of the time-varying inequality constraint are updated according to the latest window of the monitoring data to ensure that the parameters of the generated optimized intervention plan always meet the patient's current physiological state boundary.

8. The method of building a patient’s full parameter digital twin and conducting therapy simulation of claim 1, wherein, The optimized intervention parameters are input into the digital twin kernel engine for verification simulation, and a treatment simulation report is generated based on the set of verification physiological state parameters output by the simulation. The method is as follows: The optimized intervention scheme parameters are configured as the input boundary conditions or source terms of the digital twin kernel engine, driving the digital twin kernel engine to perform a verification numerical simulation based on the same physical mechanism as the original simulation but with optimized parameters, and outputting a set of verification physiological state parameters containing complete spatiotemporal evolution data. The set of verified physiological state parameters is subjected to post-processing analysis, including calculating its achievement degree with the preset clinical goal, the improvement magnitude with the predicted set of physiological state parameters, and assessing its long-term evolution trend and potential risks based on the kinetic stability criterion. The results of the post-processing analysis, together with the parameters of the optimized intervention plan and the visualization data of the key simulation process, are structurally integrated according to a preset template to generate the treatment simulation report containing quantitative assessment, plan comparison and risk warning information.

9. A system for constructing a full-parameter digital twin of a patient and performing treatment simulation, used to execute the method for constructing a full-parameter digital twin of a patient and performing treatment simulation as described in any one of claims 1-8, characterized in that, The system includes: The data integration module is configured to acquire and integrate static anatomical structure data, dynamic physiological function data, time-series metabolic parameter data, and molecular and genetic characteristic data of the target patient to construct a multi-source heterogeneous biomedical feature pool. The digital twin modeling module is configured to parse the static anatomical structure data to generate a three-dimensional anatomical structure model, and combine the dynamic physiological function data, time-series metabolic parameter data and molecular and genetic feature data to perform cross-scale fusion modeling to form the personalized digital twin initial parameter set. The multi-physics coupling simulation module is configured to drive the multi-physics coupling digital twin kernel engine with the personalized digital twin initial parameter set, simulate the effect of preset intervention measures on the patient's full-parameter digital twin, and output a set of predicted physiological state parameters. The dynamic constraint optimization module is configured to construct an inverse optimization problem based on the difference between the preset clinical goal and the predicted physiological state parameter set, and introduce dynamic physiological state constraints based on real-time monitoring data of patients during the optimization process to solve and generate optimized intervention plan parameters. The verification and report generation module is configured to input the parameters of the optimized intervention plan into the digital twin kernel engine for verification simulation, and generate a treatment simulation report based on the set of verification physiological state parameters output by the simulation.

Images