Cardiovascular blood flow limitation prediction model construction method based on hemodynamics

Abstract

The invention provides a cardiovascular blood flow limitation prediction model construction method based on hemodynamics, belongs to the technical field of cardiovascular model construction, and comprises the steps of obtaining medical image data of a to-be-tested person, synchronously collecting specific physiological parameters, reconstructing a three-dimensional geometric model and a multi-scale simulation model of a lesion blood vessel, and constructing a cardiovascular blood flow limitation prediction model. A reference model representing the blood vessel health state is generated; setting that the three-dimensional geometric model and the reference model share the same inlet boundary condition, outlet boundary condition and wall surface boundary condition; under the boundary condition, synchronous simulation calculation is conducted on the three-dimensional geometric model and the reference model through the same numerical method, and hemodynamic parameters are obtained; and respectively calculating an energy index, a functional index and a wall surface index based on the hemodynamic parameters, carrying out weighted integration on the three indexes to obtain a comprehensive score, and outputting a prediction result of blood flow limitation according to the comprehensive score.

Description

Technical Field

[0001] This invention proposes a method for constructing a cardiovascular blood flow restriction prediction model based on hemodynamics, which relates to the field of cardiovascular modeling technology. Background Technology

[0002] Cardiovascular disease is one of the leading causes of death worldwide. Its pathological basis often manifests as vascular stenosis caused by atherosclerosis, leading to myocardial or corresponding tissue ischemia. Accurately assessing the impact of vascular stenosis on hemodynamics is crucial for guiding interventional treatment decisions.

[0003] Currently, the clinical assessment of the functional significance of coronary artery stenosis mainly relies on invasive fractional flow reserve (FFR) measurement, which is limited by its invasiveness, high cost, and potential risks. Computational fluid dynamics (CFD) technology based on medical imaging can non-invasively assess hemodynamics, but existing technologies have the following shortcomings: (1) They rely solely on patient-specific models, lack quantitative comparison with healthy states, and are difficult to separate the influence of the lesion itself from other physiological factors; (2) Traditional zero-dimensional, one-dimensional, or three-dimensional single-scale models cannot take into account both the overall circulatory system and the refined analysis of local lesions; (3) The assessment indicators are singular, mostly relying only on FFR, and fail to comprehensively reflect multi-dimensional pathophysiological information such as blood flow energy dissipation and wall shear stress; (4) Boundary condition settings often use general parameters, lacking patient specificity, which affects the prediction accuracy.

[0004] Therefore, there is an urgent need for a prediction method that integrates personalized reference model comparison, multi-scale co-simulation, and multi-index comprehensive evaluation to improve the accuracy, comprehensiveness, and clinical applicability of cardiovascular blood flow restriction prediction. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of existing technologies and provide a method for constructing a cardiovascular blood flow restriction prediction model based on hemodynamics. This method constructs a patient-specific three-dimensional vascular model and its corresponding reference model, performs hemodynamic comparative simulations under unified physical conditions, and quantitatively analyzes the changes in blood flow energy caused by lesions. Simultaneously, it employs a multi-scale modeling strategy coupling zero-dimensional, one-dimensional, and three-dimensional dimensions to achieve refined simulation from macroscopic circulation to local details. Finally, it combines functional indicators and energy change indicators for comprehensive evaluation, providing accurate and comprehensive blood flow restriction prediction results for clinical practice.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: The method for constructing a cardiovascular blood flow restriction prediction model based on hemodynamics includes the following steps: Acquire medical imaging data of the subject and simultaneously collect specific physiological parameters, and reconstruct a three-dimensional geometric model of the diseased blood vessels based on the medical imaging data; A multi-scale simulation model of the three-dimensional geometric model is constructed, and a reference model representing the health status of blood vessels is generated. The reference model is obtained by performing geometric smoothing and morphological expansion operations on the narrow region in the three-dimensional geometric model to restore the diameter of the narrow region to the estimated healthy diameter. The three-dimensional geometric model is configured to share the same inlet boundary conditions, outlet boundary conditions, and wall boundary conditions with the reference model. Under the stated boundary conditions, the three-dimensional geometric model and the reference model are simultaneously simulated using the same numerical method to obtain hemodynamic parameters; Based on the hemodynamic parameters, energy, functional, and wall parameters are calculated respectively, and the three types of parameters are weighted and integrated to obtain a comprehensive score. The prediction result of blood flow restriction is output according to the comprehensive score.

[0007] Furthermore, the multi-scale simulation model is a simulation model that couples zero-dimensional, one-dimensional and three-dimensional dimensions. The zero-dimensional subsystem uses the Windkessel model to simulate the large blood vessel network and microcirculation peripheral resistance, the one-dimensional subsystem simulates pulse wave propagation based on the one-dimensional Navier-Stokes equations, and the three-dimensional subsystem constructs a fine computational domain in the lesion segment and its upstream and downstream.

[0008] Furthermore, the zero-dimensional, one-dimensional, and three-dimensional subsystems are coupled through DtN mapping and the Schwarz alternating iterative method, iteratively exchanging pressure and flow information at the interface until the convergence condition is met.

[0009] Furthermore, the estimation of the healthy diameter at the stenosis in the reference model is determined by linear interpolation or statistical extrapolation. The linear interpolation method is based on the diameter of the normal vascular segments before and after the stenosis, while the statistical extrapolation method is based on a vascular diameter-location statistical model established from a large-scale population image database.

[0010] Furthermore, the inlet boundary conditions are based on the measured flow waveform or aortic pressure waveform obtained from phase-contrast magnetic resonance imaging, and the periodic boundary conditions are reconstructed by extracting harmonics through fast Fourier transform; the outlet boundary conditions adopt the Windkessel model; and the wall boundary conditions adopt rigid wall no-slip boundary or fluid-structure interaction boundary.

[0011] Furthermore, the energy index includes energy increment and energy efficiency index, wherein the energy increment is calculated by comparing the relative percentage increase in the average total energy loss of the patient's three-dimensional geometric model and the reference model, and the energy efficiency index is calculated by the ratio of the inlet energy flow rate of the reference model to the energy loss increment of the patient model.

[0012] Furthermore, the functional indicators include fractional flow reserve, instantaneous waveformless ratio, and coronary flow reserve, which are calculated by extracting pressure time-history data and flow data from the distal end of the stenosis and the aorta in the simulation results.

[0013] Furthermore, the wall index is calculated based on the time-averaged wall shear stress, oscillatory shear index, and relative residence time. After converting each index into a dimensionless score, the highest value is taken as the wall index score.

[0014] Furthermore, the comprehensive score is calculated by weighted summation of energy indicators, functional indicators, and wall indicators, wherein the weight of each indicator is adaptively adjusted according to the target vascular bed type.

[0015] Furthermore, it also includes post-processing steps: visualizing the differences in pressure field, velocity field, and energy loss distribution between the 3D geometric model and the reference model, and automatically generating a prediction report that includes vascular geometric parameters, energy change indicators, functional indicators, comprehensive risk level, and 3D rendering.

[0016] Compared with the prior art, the present invention has the following beneficial technical effects: By constructing a reference model and comparing blood flow energy, the increase in blood flow energy dissipation caused by lesions can be accurately quantified, eliminating the interference of individual basic physiological differences and improving the specificity and sensitivity of the assessment.

[0017] Zero-dimensional, one-dimensional, and three-dimensional coupled modeling takes into account both global loops and local details, achieving an optimal balance between computational efficiency and accuracy, and is suitable for the analysis of complex vascular networks.

[0018] Combining energy indicators with traditional functional indicators provides a comprehensive reflection of hemodynamic changes, offering richer pathophysiological information and helping to identify high-risk plaques and vulnerable areas.

[0019] Non-invasive, repeatable, and low-cost, the simulation results are in high agreement with the gold standard FFR, which can provide a reliable basis for clinical decision-making in percutaneous coronary intervention (PCI) and other procedures. Attached Figure Description

[0020] To more clearly illustrate the technical solutions in the embodiments of the present invention, the drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0021] Figure 1 This is a schematic diagram of the three-dimensional geometric model of the diseased blood vessel reconstructed according to the present invention; Figure 2A schematic diagram of the coupling interface between a three-dimensional subsystem and a one-dimensional subsystem extending along the z-axis; Figure 3 This is a schematic diagram illustrating the coupling of fluid with linear elastic or hyperelastic blood vessel walls using the ALE method in this invention. Figure 4 A graph showing the relationship between time step and Courant number; Figure 5 This is a schematic diagram of the time-averaged wall shear stress (TAWSS) simulation. Detailed Implementation

[0022] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0023] In the accompanying drawings of specific embodiments of the present invention, in order to better and more clearly describe the working principle of each component in the system and show the connection relationship of each part in the device, only the relative positional relationship between each component is clearly distinguished. It does not constitute a limitation on the signal transmission direction, connection sequence, or size and shape of each part within the component or structure.

[0024] S1. Acquire medical imaging data and simultaneously collect the subject's specific physiological parameters.

[0025] S1.1 Acquisition of medical imaging data.

[0026] Cardiovascular medical imaging data are preferably obtained through the following methods: Computed tomography angiography: Dual-source CT or multi-slice spiral CT with at least 128 slices is used. Scanning parameters are set as follows: tube voltage 100-120 kVp, tube current 300-500 mAs, contrast agent injection flow rate 4-5 mL / s. The scanning phases include the early and delayed arterial phases, with an early arterial trigger threshold of 100-150 HU. The raw data are stored in DICOM 3.0 standard format, with a slice thickness not exceeding 0.5 mm and a reconstruction matrix of at least 512×512.

[0027] Magnetic resonance angiography: Three-dimensional Time-of-Flight (TOF) or contrast-enhanced magnetic resonance angiography was used, with a preferred field strength of 3.0T. A 32-channel cardiac coil was used as the receiving coil. TOF magnetic resonance angiography parameters were set to TR / TE = 25 / 3.5ms, a flip angle of 20°, and a voxel size not exceeding 0.5 × 0.5 × 0.8 mm. 3Contrast-enhanced magnetic resonance angiography used gadolinium as a contrast agent at a dose of 0.1 mmol / kg, and acquired arterial and venous phase data through dynamic tracking technology.

[0028] Digital subtraction angiography: used as the gold standard for model validation. Two-dimensional projection sequences in at least two orthogonal positions are acquired, with a frame rate of at least 15 frames per second and a pixel size not exceeding 0.2 mm.

[0029] S1.2 Simultaneously collect the following patient-specific physiological parameters.

[0030] Pressure parameters: Simultaneously record brachial artery systolic and diastolic blood pressure, and acquire real-time aortic pressure waveforms via radial or femoral artery catheters, with a sampling frequency of no less than 200Hz. For coronary arteries, distal coronary artery pressure needs to be recorded additionally.

[0031] Flow parameters: Flow waveforms in the target vessel were acquired using phase-contrast magnetic resonance imaging (PCI) or Doppler ultrasound. PCI parameters were set to a velocity encoding value of 100-150 cm / s, TR / TE = 45 / 5 ms, and a temporal resolution not exceeding 40 ms.

[0032] Blood physical properties: Venous blood samples were drawn, and blood viscosity was measured using a rotational viscometer, with a shear rate range of 0.1-1000 s⁻¹. -1 The parameters of the Carreau-Yasuda non-Newtonian model were fitted. Hemoglobin concentration, hematocrit, plasma density, and blood density were measured simultaneously, with plasma density approximately 1060 kg / m³. 3 Blood density is approximately 1065 kg / m³ 3 .

[0033] Cardiac function parameters: Left ventricular ejection fraction, heart rate, cardiac output, and RR interval on electrocardiogram are obtained by echocardiography or cardiac MRI.

[0034] S1.3 reconstructs a three-dimensional geometric model of the diseased blood vessel, such as Figure 1 As shown, it includes: Image import and noise reduction: Import DICOM data into the 3D Slicer and use Non-local Means filtering or Anisotropic Diffusion filtering to reduce image noise.

[0035] Lumen segmentation: Initial segmentation based on grayscale threshold combined with region growing algorithm is used, or U-Net / V-Net deep learning model is used for initial segmentation followed by manual correction.

[0036] Centerline extraction: The centerline of the blood vessel was extracted using the VMTK toolkit and then smoothed using B-spline processing.

[0037] 3D geometric model generation: The voxel mask is converted into a triangular patch model in STL format, and the MarchingCubes algorithm is used. Then, surface smoothing and topology repair are performed.

[0038] Model validation: The reconstructed model was registered with the digital subtraction angiography image to ensure that the geometric error RMSE does not exceed 0.3mm.

[0039] S2. Construction of multi-scale simulation models and generation of reference models for characterizing vascular health.

[0040] S2.1 Construction of multi-scale simulation model.

[0041] A zero-dimensional, one-dimensional, and three-dimensional coupled simulation model is established based on a three-dimensional geometric model. The zero-dimensional subsystem uses the Windkessel model to simulate the large vascular network and the peripheral resistance of the microcirculation, R1: R1=P m / Q mC =π•r 3 •L / 2•E•h; P m It is the average pressure, Q mC E is the average flow rate, H is the Young's modulus, h is the wall thickness, r is the radius, and L is the length.

[0042] The Windkessel model is a model used in cardiovascular hemodynamics simulation to model the pressure-flow relationship in arterial systems.

[0043] The parameters are estimated individually based on the patient's physiological data, and an inversion optimization algorithm is used to invert and optimize the one-dimensional subsystem. The pulse wave propagation is simulated based on the one-dimensional Navier-Stokes equations, using the constitutive relationship of the nonlinear elastic blood vessel wall. .

[0044] p is the pressure of a one-dimensional elastic wall. For reference area, A is the cross-sectional area of ​​the blood vessel. This represents the stiffness coefficient of the blood vessel wall. For reference pressure.

[0045] Numerical discretization employs a high-order discontinuous Galerkin method, with a CFL number ≤ 0.8. The CFL (Courant-Friedrichs-Lewy) condition is a time step constraint to ensure numerical stability. A fine computational domain is constructed for the three-dimensional subsystem within the lesion segment and its upstream and downstream sections (3-5 times the pipe diameter) using a hybrid mesh strategy. The core region mesh size is ≤ 0.1 mm, the outer region mesh size is ≤ 0.5 mm, and the initial thickness of the boundary layer expansion layer is 0.05 mm, ensuring the dimensionless wall distance parameter Y. + ≤1.0.

[0046] It should be noted that mesh generation is performed using ANSYS Meshing or OpenFOAM's snappyHexMesh.

[0047] All boundary layer meshes must meet the quality index thresholds in Table 1: ; S2.2 Generate a reference model representing vascular health Based on geometric reconstruction, stenotic regions with a diameter stenosis rate >30% are identified on a 3D geometric model. Geometric smoothing and morphological dilation operations are used to restore the diameter at the stenosis to the estimated healthy diameter. The healthy diameter is determined by interpolating the diameter of the normal vessel segment before and after stenosis or by using population statistics on diameter-length relationships, while maintaining the vessel length, branch angles, and proximal and distal anatomical structures unchanged. Based on numerical reconstruction, computational morphology methods are used to generate a smoothly transitioned vessel wall surface in the stenotic region by solving the Poisson equation, ensuring continuous curvature and no significant geometric abrupt changes. Model registration spatially registers the reference model with the 3D geometric model to ensure comparison within the same anatomical coordinate system.

[0048] It should be noted that solving the Poisson equation involves dividing the space into a regular voxel grid, typically using an octree structure to balance efficiency and accuracy. Then, methods such as the conjugate gradient method or multi-layer grid method are used to calculate the function value at each grid node. Finally, the classic moving cube algorithm is used to extract the zero isosurface from this discrete function field, forming a complete, closed, and initially smooth triangular mesh surface. Although Poisson reconstruction can generate high-quality surfaces, noise or uneven sampling in the original data can still lead to unsatisfactory local geometric details. Therefore, a geometric evolution method based on average curvature flow is further employed for fine smoothing. Each vertex on the surface moves slowly along its local normal direction, with the speed of movement determined by the curvature of the region where the point is located; the more severe the curvature, the faster the movement; flat areas remain almost stationary. This movement is mathematically equivalent to minimizing the total surface area, effectively removing high-frequency noise while preserving important large-scale morphological features.

[0049] In actual computation, for each vertex on the mesh, the algorithm examines its neighboring vertices and calculates a smoothing force based on their spatial relationships. This force essentially pulls the current vertex toward the geometric center of its local neighborhood, but is calculated using cotangent weights to ensure good numerical stability under arbitrary triangle shapes. In each iteration, all vertices synchronously move slightly along the direction of this force, forming a new surface.

[0050] To prevent excessive smoothing from flattening critical anatomical structures, the algorithm incorporates an adaptive mechanism: slowing down the movement speed in high-curvature regions and allowing for faster adjustments in flat areas. The entire iterative process continues until the average displacement of the vertices is less than a threshold related to the mesh size, or until a preset maximum number of iterations is reached, at which point the final smooth blood vessel surface is output.

[0051] It should be noted that the healthy diameter refers to the theoretical inner diameter of the target vessel segment at the same anatomical location, assuming the absence of atherosclerotic stenosis. This diameter is determined using the following two priority-based methods to ensure the anatomical accuracy of the estimated diameter: Method 1: Linear interpolation, suitable for situations where there are clearly healthy blood vessel segments at both ends of the narrow area.

[0052] In the three-dimensional geometric model, along the centerline of the blood vessel, a cross section is selected upstream of the proximal end and downstream of the distal end of the stenotic lesion. This cross section should be far away from the stenotic area and have a regular lumen shape. These sections are denoted as the proximal reference cross section and the distal reference cross section, respectively.

[0053] Measure the equivalent diameter D of these two sections. prox With D dist .

[0054] For any point P on the centerline within the narrow region, let the arc length of the centerline from P to the center point of the near-end reference section be L. p The arc length of the centerline to the center point of the far-end reference section is L. d The total arc length between the two points is L. total Then the healthy diameter D at point P is... healthy (P) is determined by linear interpolation: ; This method is based on the physiological assumption that the diameter of blood vessels usually changes gradually along the centerline, which can preserve the individual vascular anatomy of the patient to the greatest extent.

[0055] It should be noted that linear interpolation is only applicable when the following conditions are met simultaneously: A. Distance requirements: The centerline distance between the near-end reference section and the near-end boundary of the narrow region is ≥5mm, and the centerline distance between the far-end reference section and the far-end boundary of the narrow region is ≥5mm.

[0056] B. Shape Requirements: The reference section should meet the following requirements: The standard deviation of the lumen diameter is ≤ 10% of the average diameter; Lumen roundness ≥ 0.85: There was no obvious local dilation or stenosis.

[0057] C. Continuity requirement: The vascular segments at both ends of the narrow area should maintain anatomical continuity, without obvious branching or abrupt changes in anatomical structure.

[0058] If any of the following conditions are not met, method two should be used: 1. The distance between the reference section and the narrow area is <5mm; 2. The roundness of the reference section is <0.85; 3. The standard deviation of the reference section diameter is >10% of the average diameter; 4. The diameter variation rate at both ends of the narrow area is ≥30%; 5. The length of the narrow area is >10 times the diameter of the blood vessel.

[0059] Method 2: Statistical extrapolation, which is suitable for narrow, diffuse, or difficult-to-determine reference sections.

[0060] When Method 1 is not applicable, a statistical model of blood vessel diameter-location based on a large-scale population image database is used for estimation.

[0061] Specifically, statistical extrapolation is applicable to the following situations: A: The length of the stenotic region is ≥ 10 times the diameter of the blood vessel; B: The diameter change rate at both ends of the narrow region is ≥30%; C: The narrow area contains multiple narrow points; D: Unable to obtain a reference section that meets the distance and shape requirements.

[0062] The statistical model is derived from a publicly available or self-constructed database covering healthy individuals. Database requirements: The statistical model must be built based on CTA or MRA data from at least 100 healthy blood vessels, covering the target blood vessel type, such as LAD, LCX, RCA, etc.

[0063] By analyzing the correlation between the diameter of a specific blood vessel and the arc length from the centerline of the coronary artery ostium, a univariate linear regression model was established: D expected =a×L+b, where L is the arc length from the opening, and a and b are the population specificity coefficients obtained through regression analysis.

[0064] By substituting the arc length L from the coronary artery ostium at each point on the centerline of the target lesion segment into the statistical model above, the expected healthy diameter D corresponding to each point can be calculated. expected .

[0065] In obtaining the healthy diameter D expected After spatial distribution, using this as the target, the diameter of the narrow section is restored to the estimated healthy diameter through geometric smoothing and morphological expansion operations.

[0066] Preferably, the surface of the vessel wall in the target area is smoothed using an iterative smoothing algorithm based on curvature flow. In each iteration, the average curvature of the adjacent faces of each vertex is calculated, and the vertex is moved along the normal to a position with more uniform curvature. The number of iterations is dynamically adjusted according to the length of the target area, with three to five iterations corresponding to each millimeter of target area length, and the total number of iterations is controlled within fifty. The process automatically stops when the maximum curvature change of the surface is less than five percent after two consecutive iterations. During the smoothing process, the total length of the vessel, the branch angle, and the proximal and distal anatomical structures are kept strictly unchanged, and constraints are achieved by fixing the coordinates of non-target area vertices.

[0067] In the smoothed target area, the lumen is expanded layer by layer along the normal direction of the vessel centerline. The expansion amplitude is adaptively adjusted point by point to make the diameter at each point approximate the estimated healthy diameter value. The expansion step size is set to 10% to 20% of the difference between the current diameter and the healthy diameter to avoid geometrical abrupt changes caused by a single expansion. A buffer zone is used at the boundary between the target area edge and the normal vessel segment for a smooth transition. The width of the buffer zone is 0.5 mm extending from the target area edge to both sides. Within this area, the expansion amplitude linearly decreases to zero to ensure continuous curvature and no step-like geometric defects.

[0068] The generated reference model was spatially registered with the original patient model, and an iterative nearest-point algorithm was used to ensure that the centerlines of the two models overlapped by more than 99%. After registration, the diameter differences between the two models were compared at the same anatomical location to ensure that the deviation between the estimated target area diameter of the reference model and the healthy diameter did not exceed 0.1 mm. Finally, the surface of the reference model was free of self-intersections, and the aspect ratio of the triangular facets was controlled within 5:1.

[0069] S2.3 Multiscale Coupling The three-dimensional and one-dimensional coupling is achieved by using DtN mapping at the inlet and outlet sections of the three-dimensional geometric model. The convergence tolerance c of pressure and flow rate is exchanged through the Schwarz alternating iteration method and is a dimensionless quantity.

[0070] Specifically, for microvessels with a diameter <2.0 mm, due to the weak flow signal in small vessels, excessively strict tolerances can easily lead to divergence; therefore, the convergence tolerance c can be relaxed to 1.5 × 10⁻⁶. -3 ; For blood vessels with a diameter >3.5 mm, the convergence tolerance c can be tightened to 8.0 × 10⁻⁶ because larger vessels are more sensitive to changes in flow rate. -4 .

[0071] One-dimensional and zero-dimensional coupling is performed on a one-dimensional model, and the terminal uses the feature line method to match the zero-dimensional Windkessel model.

[0072] Specifically, such as Figure 2 As shown, in the three-dimensional subsystem With a one-dimensional subsystem extending along the z-axis The coupling interface, namely the inlet and outlet sections of the three-dimensional geometric model, hereinafter referred to as the interface Γ: the three-dimensional geometric model provides velocity, and the one-dimensional model provides pressure, which are exchanged through DtN (Dirichlet-to-Neumann) mapping; at the coupling interface between the one-dimensional subsystem and the zero-dimensional subsystem, namely the one-dimensional model terminal: the one-dimensional model provides flow and pressure characteristics, and the zero-dimensional model provides impedance boundary conditions.

[0073] Preferably, a strongly coupled solution between subsystems is achieved through DtN mapping and the Schwarz alternating iterative method: An interface Γ is established at the inlet and outlet sections of the 3D model, and bidirectional pressure-flow rate transfer is achieved through iterative exchange. Specifically, at the interface Γ, iteration is performed at each time step: First, the current estimated flow rate is used as the Neumann condition for the 3D model to obtain the interface pressure distribution, which is then used as the Dirichlet condition for the 1D model after area-weighted averaging; subsequently, the 1D model is solved to obtain the new interface flow rate. This process is repeated until the relative changes in interface flow rate and pressure are both less than the convergence tolerance c or the iteration upper limit N is reached. This process realizes bidirectional real-time exchange of pressure and flow rate between the 3D and 1D models.

[0074] The upper limit N for the number of iterations is set based on the following: tests on 500 clinical vascular models, including 10 typical stenosis morphologies, showed that 97.3% of the cases converged within 60 iterations, with a convergence rate of 97.3%, and 99.8% of the cases converged within 100 iterations. The convergence rates are shown in Table 2.

[0075] Table 2 ; In a preferred embodiment, when the iteration count reaches N and still has not converged, the input geometric data is checked to verify whether the radius of curvature of the vessel centerline is greater than or equal to the vessel diameter. If the radius of curvature is less than the diameter, it is marked as geometric discontinuity. The boundary conditions are checked to confirm that the inlet flow rate matches the pressure. If the inlet flow rate is greater than 120% of the outlet flow rate, it is marked as boundary abnormality.

[0076] Hierarchical processing: First-level processing: Temporarily increase N to N+N / 2, and rerun the iteration.

[0077] Secondary processing: If convergence is still not achieved, the result of the previous iteration is used as the final solution, and a convergence warning is marked in the output file.

[0078] Level 3 processing: If two consecutive retries fail, the simulation is paused and a diagnostic report is output, requiring the user to manually correct and resubmit.

[0079] At the terminals of the one-dimensional model and the zero-dimensional model, the method of characteristics is used for coupling: by using the forward characteristic variables of the one-dimensional model terminal and the terminal pressure calculated by the zero-dimensional model, the matching terminal flow rate is obtained by solving the characteristic line compatibility relationship, thereby realizing the non-reflective connection between the one-dimensional model and the downstream circulating impedance.

[0080] The S3-constrained 3D geometric model and the reference model share the same entry conditions.

[0081] S3.1 Inlet boundary conditions.

[0082] The three-dimensional geometric model and the reference model share the same entry conditions. The measured flow waveform Q(t) from PC-MRI or the aortic pressure waveform Pa(t) is used. The periodic boundary turbulence intensity is set to 5-10% by extracting the first 10 harmonics through FFT.

[0083] S3.2 Export Boundary Conditions The zero-dimensional model outlet is connected to the Windkessel model, with the terminal venous pressure set to 5 mmHg. The three-dimensional geometric model outlet uses a mass flow rate outlet at the three-dimensional-one-dimensional coupling interface, with the flow rate value provided in real time by the one-dimensional model.

[0084] S3.3 Wall Boundary Conditions Rigid walls are assumed to be non-slip walls (u=0) by default. Fluid-structure interaction (ALE) is optional, using the fluid-structure interaction method to couple the fluid to the linearly elastic or hyperelastic vessel wall. Figure 3 As shown.

[0085] S4. Simulate the two models simultaneously under the same numerical methods and boundary conditions.

[0086] Regarding the solver settings in S4.1, the PISO or SIMPLEC algorithm is used to solve the unsteady Navier-Stokes equations. Second-order Backward Euler scheme is used for time discretization, with a time step Δt ranging from 0.001 to 0.01 seconds to ensure the Courant number does not exceed 0.5. Figure 4 As shown. The choice of blood model depends on the shear rate: when the Newtonian fluid assumption is met, a Newtonian fluid model with a viscosity of 0.0035 Pa·s is used; if the shear rate is low or the shear thinning effect needs to be considered, the Carreau-Yasuda non-Newtonian model is selected to more accurately describe the blood rheological properties.

[0087] Define the threshold for shear rate: When the shear rate γ is greater than the shear rate threshold, such as 100s -1 When choosing, select the Newtonian fluid model.

[0088] When the shear rate γ is less than or equal to the shear rate threshold, a non-Newtonian fluid model, such as the Carreau-Yasuda model, is selected.

[0089] For Newtonian fluid models, viscosity μ can be obtained experimentally, for example, by using a rotational viscometer to measure the viscosity of blood at high shear rates.

[0090] For non-Newtonian fluid models, parameters such as zero-shear viscosity, infinite-shear viscosity, and relaxation time can be obtained by fitting experimental data.

[0091] Example of blood viscosity fitting: Patient data: Male, 45 years old, fasting blood glucose 6.2 mmol / L, measurement data are shown in Table 3: Table 3 ; In the S4.2 comparative simulation execution phase, efficient parallel computing is first achieved through domain decomposition and MPI parallel technology, and GPU acceleration is used to solve the Poisson pressure equation to address computational bottlenecks. Next, synchronous initialization is performed, unifying the initial pressure across the entire field of both comparative models to the mean diastolic pressure, while simultaneously initializing the velocity field of all nodes to zero. To achieve a physically realistic flow field state, steady-state startup calculations are performed for 3 to 5 cardiac cycles to allow the fluid to fully develop and reach a periodic steady state. After startup, the transient comparative calculation phase begins, running both models synchronously for at least 10 cardiac cycles. Data acquisition starts from the 6th cycle, recording detailed flow field data for the last 5 complete cycles for subsequent analysis.

[0092] S4.3 Blood flow energy calculation.

[0093] Blood flow energy consists of three parts: pressure potential energy, kinetic energy, and viscous dissipation energy. At the inlet section of the model, the instantaneous energy flow rate is calculated using surface integrals, and the integral calculation must include pressure and velocity-related terms. At all outlet sections, the same integral method is used to calculate the instantaneous energy flow rate at each outlet. Finally, energy loss is determined by subtracting the sum of the energy flow rates at all outlets from the inlet energy flow rate, and the energy loss over the entire cardiac cycle is time-averaged to obtain a quantitative index of average energy dissipation over that cycle.

[0094] In practical applications, the calculation methods for instantaneous energy flow rate and energy loss are as follows: Based on the instantaneous pressure and velocity field data obtained from hemodynamic simulations, the energy flow rate is calculated using an area integral method across any cross-section of the blood vessel. Specifically, the product of the pressure value and the flow velocity at each point on the cross-section is integrated to obtain the pressure potential energy term; simultaneously, the cube values ​​of the blood density and flow velocity are integrated and multiplied by half to obtain the kinetic energy term; the sum of these two is the instantaneous energy flow rate of that cross-section. Performing the above calculations on the inlet cross-section and all outlet cross-sections of the three-dimensional geometric model yields the inlet and outlet instantaneous energy flow rates. Subtracting the sum of all outlet energy flow rates from the inlet energy flow rate gives the instantaneous energy loss. Based on this, the time-averaged instantaneous energy loss is calculated for each moment within a complete cardiac cycle, ultimately yielding the time-averaged total energy loss, expressed in watts.

[0095] S5. Calculate energy, functional, and wall surface indicators, integrate them according to weights to calculate a comprehensive score, form a clinical grade, and generate a final report including a 3D rendering.

[0096] S5.1, Energy Indicators Energy indicators include energy increment and energy efficiency index.

[0097] The energy increment is calculated by comparing the average total energy loss of the patient's three-dimensional geometric model with that of the reference model. The energy efficiency index is calculated by the ratio of the average energy flow rate at the inlet of the reference model to the energy loss increment of the patient's three-dimensional geometric model. The pressure drop energy ratio is the ratio of pure pressure potential energy loss to total energy loss, used to quantify the contribution of pressure drop to total energy loss.

[0098] Energy increment is defined as the relative percentage increase between the average total energy loss of the patient-specific three-dimensional geometric model and the average total energy loss of the reference model. First, the average total energy loss of the patient model and the reference model within the same cardiac cycle is calculated separately. Then, the difference in average energy loss between the two models is divided by the average energy loss of the reference model, and the result is multiplied by 100%. When the calculated energy increment exceeds 50%, it is considered a high-risk energy dissipation state, indicating severe blood flow restriction.

[0099] The Energy Efficiency Index (EEI) is defined as a correction term for the ratio between the inlet energy flow rate of the reference model and the energy loss increment of the patient model. It is used to quantify the degree of decline in blood flow energy transfer efficiency caused by lesions. Specifically, the average time-to-time energy flow rate of the reference model at the inlet is first determined. Then, the difference between the average time-to-time energy loss of the patient model and the average time-to-time energy loss of the reference model is calculated, i.e., the energy loss increment. This energy loss increment is divided by the inlet energy flow rate of the reference model, and then the quotient is subtracted by one to obtain the EEI. The EEI value ranges from zero to one. When the EEI is less than 0.7, it indicates that the lesion has led to a significant reduction in energy transfer efficiency, and the risk of blood flow restriction is high.

[0100] The two indicators mentioned above are strictly based on the law of conservation of energy. The energy increment directly reflects the additional mechanical work of blood flow consumed by the diseased vessel relative to the healthy state. This increment mainly originates from enhanced viscous dissipation and turbulent losses in the narrowed region. The energy efficiency index comprehensively considers the impact of the lesion on the overall energy transmission efficiency. The threshold settings for these two indicators are based on the following: receiver operating characteristic (ROC) curve analysis using retrospective computational fluid dynamics simulation data and the gold standard measurement of invasive fractional flow reserve from one hundred patients diagnosed with coronary artery disease. When the energy efficiency index is less than 0.7, its diagnostic specificity reaches 85%, and its sensitivity reaches 82%, providing clear statistical evidence.

[0101] Energy increment and energy efficiency index quantify the impact of lesions from the perspective of blood flow mechanical work loss, but a single energy indicator cannot comprehensively reflect the microcirculatory status, myocardial oxygen supply matching, and pathophysiological changes of the vascular wall. Therefore, it is necessary to simultaneously extract traditional functional pressure-flow indicators and wall shear stress indicators for multi-dimensional cross-validation.

[0102] S5.2, Functional Indicators.

[0103] The simulation results extract the time-series data of pressure distal to the stenosis and aortic pressure. The ratio of the average pressure of the two under the congested state is used to obtain the fractional flow reserve. When the fraction is less than or equal to 0.80, the blood flow function is considered restricted. The ratio of the average pressure of the two during the diastolic waveless period under the resting state is used to obtain the instantaneous waveless ratio. When the ratio is less than or equal to 0.89, the blood flow reserve is considered restricted. The ratio of peak flow in the congested state to peak flow in the resting state is used to obtain the coronary flow reserve (CFR). When the ratio is less than 2.0, it indicates abnormal microcirculation function.

[0104] S5.3, Wall surface index.

[0105] TAWSS, OSI, and RRT need to be converted into dimensionless wall metrics. The specific construction method is as follows: For the time-averaged wall shear stress (TAWSS), identify regions in all mesh elements on the lumen surface where the TAWSS is less than 0.4 Pa, and calculate the percentage of this region's area relative to the total lumen area. Figure 5 As shown. When the percentage is zero, the TAWSS sub-item score is zero; when the percentage is 10% or more, the TAWSS sub-item score is 100; when it is in between, a linear score is given by multiplying the actual percentage by 100.

[0106] For the oscillating shear index (OSI), grid cells with an OSI greater than 0.2 are identified, and their total area is calculated as a percentage of the total lumen area. The scoring rules are the same as for the TAWSS sub-item: zero points for a percentage of 0, one hundred points for a percentage of 10% or higher, and linear interpolation for intermediate states.

[0107] For the relative residence time (RRT), calculate the logarithmic mean of the RRT values ​​of all grid cells. A score of zero is given when the logarithmic mean is less than one, and a score of one hundred is given when it is greater than or equal to five. For values ​​between one and five, a linear transformation is performed proportionally.

[0108] After obtaining the scores for the three sub-items, the wall index is selected based on the highest value among the three. This design is based on the conclusions of clinical fluid dynamics research: if any of the indicators in TAWSS, OSI, and RRT reaches the high-risk threshold, it indicates a significant decrease in plaque stability. Therefore, the maximum value principle is used instead of a simple average to avoid underestimating the risk.

[0109] S5.3 Comprehensive Assessment and Grading.

[0110] The three categories of indicators are weighted and summed to obtain the comprehensive score S:

[0111] in, , , The scores are for the normalized energy index, functional index, and wall index, respectively.

[0112] Energy index weight The preferred percentage is 30%. A high risk is defined as ΔE% greater than 50% or EEI less than 0.7%.

[0113] Functional indicator weights The preferred value is 50%: the main criterion is that FFR is less than or equal to 0.80 or iFR is less than or equal to 0.89, where iFR is the instantaneous waveform-free ratio.

[0114] Wall surface index weight The preferred percentage is 20%: when the area of ​​a region with an OSI greater than 0.2 and a TAWSS less than 0.4 is greater than 10%, it is considered a high-risk plaque marker.

[0115] The weighting of energy, function, and wall parameters is not fixed but adaptively adjusted based on the pathophysiological characteristics of the target vascular bed, hemodynamic response patterns, and clinical risk assessment priorities. The adjustment mechanism is based on a subgroup analysis of a multicenter clinical dataset including coronary arteries, carotid arteries, renal arteries, and lower extremity arteries. This dataset encompasses 1,200 cases of blood flow restriction diagnosed by imaging and functional examinations. Specific adjustment rules are as follows: Coronary artery weighting: This is the default setting used for most coronary artery lesion assessments. Functional indicators are maintained at the highest priority (50%) because myocardial oxygen supply directly depends on coronary pressure-flow matching, and the predictive value of FFR and iFR for myocardial ischemia events has been confirmed by large-scale randomized controlled trials. Energy indicators account for 30% of the weighting, used to capture the additional mechanical energy consumption caused by stenosis, and have supplementary diagnostic value for diffuse and tandem lesions. Wall indicators account for 20% of the weighting, primarily identifying vulnerable plaques and providing auxiliary information for risk stratification of acute coronary syndromes.

[0116] Carotid artery weighting: When applied to carotid artery stenosis assessment, the weight of wall parameters is increased to 35%, the weight of functional parameters is correspondingly decreased to 40%, and the weight of energy parameters remains at 25%. This adjustment is based on the fact that the primary clinical risk of carotid artery stenosis stems from thromboembolic events caused by plaque rupture, rather than simply insufficient cerebral perfusion. Therefore, the oscillatory shear index and time-averaged wall shear stress are more important for assessing plaque fibrous cap stability. Although the weight of functional parameters is reduced, they remain core indicators because the carotid pressure gradient directly affects cerebral blood flow reserve.

[0117] Renal artery weighting: When assessing renal artery stenosis, the weight of energy indicators is increased to 40%, the weight of functional indicators is decreased to 45%, and the weight of wall indicators is 15%. This adjustment is because the kidneys are extremely sensitive to decreased perfusion pressure; energy dissipation caused by renal artery stenosis directly reflects glomerular filtration pressure loss, which is associated with deteriorating renal function and refractory hypertension. Functional indicators remain the primary weighting, but the pressure ratio threshold needs to be adjusted to a renal artery-specific standard. The weight of wall indicators is reduced because renal artery atherosclerotic plaques are relatively difficult to rupture and therefore receive less clinical attention.

[0118] Implementation of weight adjustment: Before clinical application, the system automatically calls the corresponding weight configuration based on the target vascular bed type input by the user. The weight adjustment range is determined based on the meta-analysis results, ensuring that the area under the receiver operating characteristic curve (AUC) of the adjusted composite score is not less than 0.9. If the assessment subject has mixed vascular lesions or multiple sites involved, the system supports manual fine-tuning of the weights. The fine-tuning range is limited to within ±20% of the default value. If the range is exceeded, the report must indicate that it is a non-standard configuration and warn that the interpretation of the results should be approached with caution.

[0119] Comprehensive rating: Unrestricted: All indicators are normal.

[0120] Mildly restricted: FFR is between 0.80 and 0.85 and includes 0.85, while ΔE% is between 20% and 40% and includes 40%.

[0121] Moderately restricted: FFR is between 0.75 and 0.80 and includes 0.80, or iFR is between 0.85 and 0.89 and includes 0.89, while ΔE% is between 40% and 60% and includes 60%.

[0122] Severe limitation: FFR less than or equal to 0.75, or ΔE% greater than 60%, or EEI less than 0.6%.

[0123] S5.4 Post-processing and Reporting Visualization: The pressure field and velocity field differences between the 3D geometric model and the reference model are overlaid and displayed in ParaView software, and the energy loss distribution is presented using a pseudo-color cloud map.

[0124] Report generation: Automatically generates PDF reports containing vascular geometric parameters, energy change indicators, functional indicators, comprehensive risk level, and 3D renderings.

[0125] Example 2 Take a 65-year-old male patient with 70% proximal LAD stenosis as an example: Three-dimensional geometric model: CTA reconstruction, stenosis rate 68%, lesion length 15mm.

[0126] Reference model: The diameter of the narrow section was restored to 3.5 mm (originally 1.2 mm) through geometric smoothing.

[0127] Simulation: The two models ran synchronously for 10 cycles, and the calculation took about 50 minutes.

[0128] Results: The time-averaged energy loss of the 3D geometric model is E=0.85W, while that of the reference model is E=0.32W, with ΔE%=166%; FFR=0.73, iFR=0.84, and the area with OSI>0.2 accounts for 18%. Overall, the system is considered severely limited, and immediate PCI is recommended. The measured invasive FFR is 0.74, with a simulation error of 1.3%, indicating good consistency between energy and performance indicators.

[0129] ; In vitro experiments: 3D printing of personalized and reference models, PIV measurement of velocity field, comparison with CFD results, and correlation coefficient R of velocity distribution. 2 ≥0.85.

[0130] Clinical trial: A prospective study of 100 patients was conducted. The difference between the simulated FFR and the invasive FFR was <0.02, with a 95% consensus margin of ±0.05 and an area under the ROC curve (AUC) ≥0.92. The energy change index ΔE% was significantly negatively correlated with FFR (r=-0.78, p<0.001).

[0131] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope disclosed in this application, and these modifications or substitutions should all be covered within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for constructing a cardiovascular blood flow restriction prediction model based on hemodynamics, characterized in that, Includes the following steps: Acquire medical imaging data of the subject and simultaneously collect specific physiological parameters, and reconstruct a three-dimensional geometric model of the diseased blood vessels based on the medical imaging data; A multi-scale simulation model of the three-dimensional geometric model is constructed, and a reference model representing the health status of blood vessels is generated. The reference model is obtained by performing geometric smoothing and morphological expansion operations on the narrow region in the three-dimensional geometric model to restore the diameter of the narrow region to the estimated healthy diameter. The three-dimensional geometric model is configured to share the same inlet boundary conditions, outlet boundary conditions, and wall boundary conditions with the reference model. Under the stated boundary conditions, the three-dimensional geometric model and the reference model are simultaneously simulated using the same numerical method to obtain hemodynamic parameters; Based on the hemodynamic parameters, energy, functional, and wall parameters are calculated respectively, and the three types of parameters are weighted and integrated to obtain a comprehensive score. The prediction result of blood flow restriction is output according to the comprehensive score.

2. The method according to claim 1, characterized in that, The multi-scale simulation model is a coupled simulation model of zero-dimensional, one-dimensional and three-dimensional components. The zero-dimensional subsystem uses the Windkessel model to simulate the large vascular network and microcirculation peripheral resistance. The one-dimensional subsystem simulates pulse wave propagation based on the one-dimensional Navier-Stokes equations. The three-dimensional subsystem constructs a fine computational domain in the lesion segment and its upstream and downstream.

3. The method according to claim 2, characterized in that, The zero-dimensional, one-dimensional, and three-dimensional subsystems are coupled through DtN mapping and the Schwarz alternating iterative method, and pressure and flow information are iteratively exchanged at the interface until the convergence condition is met.

4. The method according to claim 1, characterized in that, The healthy diameter at the stenosis in the reference model is estimated using either linear interpolation or statistical extrapolation. Linear interpolation is based on the diameter of the normal vascular segments before and after the stenosis, while statistical extrapolation is based on a vascular diameter-location statistical model established from a large-scale population image database.

5. The method according to claim 1, characterized in that, The inlet boundary conditions are based on the measured flow waveform or aortic pressure waveform obtained from phase-contrast magnetic resonance imaging, and the periodic boundary conditions are reconstructed by extracting harmonics through fast Fourier transform; the outlet boundary conditions are based on the Windkessel model; and the wall boundary conditions are based on rigid wall no-slip boundary or fluid-structure interaction boundary.

6. The method according to claim 1, characterized in that, The energy indicators include energy increment and energy efficiency index. The energy increment is calculated by comparing the relative percentage increase in the average total energy loss of the patient's three-dimensional geometric model and the reference model. The energy efficiency index is calculated by the ratio of the inlet energy flow rate of the reference model to the energy loss increment of the patient model.

7. The method according to claim 1, characterized in that, The functional indicators include fractional flow reserve, instantaneous waveformless ratio, and coronary flow reserve, which are calculated by extracting pressure time-history data and flow data from the distal end of the stenosis and the aorta in the simulation results.

8. The method according to claim 1, characterized in that, The wall index is calculated based on the time-averaged wall shear stress, oscillatory shear index, and relative residence time. After converting each index into a dimensionless score, the highest value is taken as the wall index score.

9. The method according to claim 1, characterized in that, The comprehensive score is calculated by weighted summation of energy indicators, functional indicators, and wall indicators, with the weight of each indicator adaptively adjusted according to the target vascular bed type.

10. The method according to claim 1, characterized in that, It also includes post-processing steps: visualizing the differences in pressure field, velocity field, and energy loss distribution between the 3D geometric model and the reference model, and automatically generating a prediction report containing vascular geometric parameters, energy change indicators, functional indicators, comprehensive risk level, and 3D rendering.

Images