Simulation system for evaluating cardiac ablation plans

Abstract

The present application provides a simulation system for evaluating cardiac ablation schemes, including: an image acquisition module for acquiring multiple chest images; a three-dimensional grid reconstruction module for marking ablation points on the multiple chest images and obtaining a three-dimensional grid of the heart by a finite element method; a fiber rotation direction calculation module for calculating the myocardial fiber rotation direction based on the three-dimensional grid of the heart, and obtaining the direction vector of the myocardial fiber rotation direction corresponding to each grid point on the three-dimensional grid of the heart; an electrophysiological parameter simulation module for simulating the cardiac action potential excitation process based on electrophysiological parameters, the geometric positions of grid points, and the fiber rotation directions corresponding to the grid points; a dynamic simulation module for simulating cardiac electrical signals based on the cardiac action potential excitation process to obtain dynamic simulation results. The present application introduces anisotropy through fiber rotation direction, making the simulation results more realistic.

Description

Technical Field

[0001] The present application mainly relates to the field of biomedical engineering simulation, and in particular to a simulation system for evaluating cardiac ablation schemes. Background Art

[0002] Some types of arrhythmias are caused by local abnormal areas of the heart's electrical system. Radiofrequency ablation can usually be used to intervene in the abnormal area to eliminate the arrhythmia. In order to determine the ablation point, it is usually necessary to use an invasive ECG detection device for pacing stimulation, detect electrophysiological signals, and then determine the ablation point based on the electrophysiological signals. At present, due to the limitations of technology and principles, the effect of ablation is often uncertain, and arrhythmias have a large recurrence rate. For recurrent patients, the aforementioned invasive pacing and detection need to be performed again, which is undoubtedly very painful for patients. At present, there are also some technical paths for the development of non-invasive electrophysiology. Usually, more electrodes are used to collect multi-lead ECG signals on the surface of the trunk, and then the distribution of the epicardial action potential phase diagram is inferred in the form of an ECG inverse problem through the conduction matrix, and then the lesions and ablation plans are evaluated.

[0003] However, the related technologies of the ECG inverse problem face difficulties in principle and accuracy. In principle, the ECG inverse problem relies on the boundary element method to calculate the conduction matrix between the geometric grids of the epicardium and the trunk surface to establish the relationship between the electric potentials on the corresponding grids. The calculation of the conduction matrix is ​​extremely sensitive to the accuracy of the geometric grid and is interfered by the physiological parameters of the patient's own chest, lungs and other organs, resulting in limited accuracy of the obtained conduction matrix. In addition, the equations solved by the inverse problem are severely underdetermined, which makes the numerical solution seriously unstable during calculation. The accuracy and number of electrodes on the trunk surface as input parameters can never reach the accuracy required to eliminate this instability for the underdetermined equation.

[0004] Therefore, the accuracy of the solution obtained by the inverse problem is very limited. If it is used to predict the arrhythmia lesion, its accuracy is not reliable, so this non-invasive assessment of the ablation point cannot accurately obtain the location of the lesion. Summary of the invention

[0005] In response to the above technical problems, the present application provides a simulation system that can more accurately simulate the actual electrophysiological activity of myocardial fibers.

[0006] In order to solve the above technical problems, the present application provides a simulation system for evaluating cardiac ablation schemes, including: an image acquisition module, used to acquire multiple chest images, wherein the chest images include cardiac angiography images; a three-dimensional grid reconstruction module, used to mark ablation points on the multiple chest images, and obtain a three-dimensional grid of the heart by a finite element method, wherein the ablation points include estimated ablation points to be ablated, and / or actual ablation points where ablation has been performed; a fiber rotation direction calculation module, used to calculate the myocardial fiber rotation direction based on the three-dimensional grid of the heart, and obtain the myocardial fiber rotation direction corresponding to each grid point on the three-dimensional grid of the heart. A direction vector of fiber rotation; an electrophysiological parameter simulation module, which is used to simulate the cardiac action potential excitation process according to the electrophysiological parameters, the geometric positions of the grid points, and the fiber rotation directions corresponding to the grid points, wherein the electrophysiological parameters include cell-level parameters and tissue-level parameters, the cell-level parameters include current parameters of ion channels, and the tissue-level parameters include the conductivity of different heart regions, pacing stimulation parameters, and the propagation rate of electrical signals set according to the direction vector of the myocardial fiber rotation direction; a dynamic simulation module, which is used to simulate the cardiac electrical signals according to the cardiac action potential excitation process to obtain dynamic simulation results.

[0007] In one embodiment of the present application, the chest image includes an LGE-MRI enhanced image.

[0008] In one embodiment of the present application, the three-dimensional grid reconstruction module is further used to segment the multiple chest images, and the segmented boundaries include the left ventricular endocardium, the right ventricular endocardium and the epicardium.

[0009] In one embodiment of the present application, the fiber handedness calculation module calculates the myocardial fiber handedness according to the three-dimensional heart mesh based on the LD-RB algorithm.

[0010] In one embodiment of the present application, the electrophysiological parameter simulation module obtains the current parameters of the ion channel corresponding to each grid point on the three-dimensional grid of the heart based on the ion channel model.

[0011] In one embodiment of the present application, the electrophysiological parameter simulation module determines the conductivity of different regions in the three-dimensional grid of the heart based on the segmented chest image, wherein the conductivity corresponding to the ablation point is different from the conductivity corresponding to normal myocardial cells.

[0012] In one embodiment of the present application, the conductivity corresponding to the ablation point is 0.

[0013] In one embodiment of the present application, the pacing stimulation parameters include parameters corresponding to a pre-programmed stimulation scheme, and the pre-programmed stimulation scheme is configured to include:

[0014] Step 1: Select the location of the stimulation point;

[0015] Step 2: Use S1S2 stimulation, the stimulation sequence includes 6 consecutive stimulations of the stimulation point by the pacing stimulation S1 with a cycle length of T1, and after the last S1 stimulation, apply an S2 stimulation, the cycle length of the S2 stimulation is T2, if no abnormality is induced after the S2 stimulation, then the cycle length T2 of the S2 stimulation is decreased according to a preset step length, and the S2 stimulation is applied cyclically until the abnormality is induced;

[0016] Step 3: If the S2 stimulation still fails to induce abnormality, then the S3 stimulation is applied, and the cycle length T3 of the S2 stimulation is decreased according to the preset step length, and the S3 stimulation is applied cyclically until the abnormality is induced;

[0017] Step 4: If the S3 stimulation still fails to induce abnormality, then the S4 stimulation is applied, and the cycle length T4 of the S4 stimulation is decreased according to the preset step length, and the S4 stimulation is applied cyclically until the abnormality is induced;

[0018] Step 5: If S4 stimulation still fails to induce abnormalities, it is considered that the stimulation point cannot induce the desired arrhythmia.

[0019] In one embodiment of the present application, the dynamic simulation module simulates the cardiac electrical signal based on the Reaction-Diffusion diffusion process, and the simulation process includes:

[0020] Step 10: The transmembrane position and reset voltage on the time-space evolution scale are used as contraction variables u, v and nonlinear response functions f, g to satisfy the equation:

[0021]

[0022] Among them, D1 and D2 are the coefficients corresponding to the shrinkage variables u and v, is the Laplace operator,

[0023] Step 20: Discretize the equation in step 10 in the time dimension, set t = 1, 2, ..., T, and record the grid point sequence i = 1, 2, ..., N in the three-dimensional grid of the heart, and obtain the following equation:

[0024]

[0025] Step 30: Write the equation in step 20 as the following matrix:

[0026]

[0027] in,

[0028] Step 40: Use the optimization method to solve the matrix in step 30 and obtain U t and V t .

[0029] In one embodiment of the present application, the Fitzhugh-Nagumo model is used in step 10, so that:

[0030] f(u,v)=C1u(1-u)(ua)-C2uv

[0031] g(u,v)=b(u-dv)

[0032] Among them, C1 and C2 are parameters related to the curve shape of the function f, 1 / b is the average time constant of the slow current passing through the ion channel, a is a parameter related to the resting voltage, and d is the ratio of the transmembrane current to the slow local recovery current in the ion channel.

[0033] In one embodiment of the present application, it further includes: an evaluation module, which is used to evaluate the cardiac ablation plan according to the dynamic simulation result, and the evaluation result includes whether to change the ablation point.

[0034] The simulation system of the present application constructs a three-dimensional cardiac grid through images, calculates the rotation direction of myocardial fibers, obtains the cardiac action potential excitation process with fiber rotation anisotropy, introduces the anisotropic cardiac potential diffusion process in the dynamic simulation module for simulation calculation, simulates the actual electrophysiological properties of myocardial fibers, greatly improves the numerical accuracy, and avoids the instability of the results obtained by the boundary element method, circumvents the defects of the ECG inverse problem, and can simulate the actual case characteristics more finely. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] The accompanying drawings are included to provide a further understanding of the present application. They are included and constitute a part of the present application. The accompanying drawings illustrate embodiments of the present application and together with the present specification serve to explain the principles of the present application. In the accompanying drawings:

[0036] Figure 1 is a block diagram of a simulation system for evaluating cardiac ablation schemes of the present application;

[0037] Figure 2 This is an exemplary flow chart of using the simulation system according to an embodiment of the present application for evaluating a cardiac ablation plan. DETAILED DESCRIPTION

[0038] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the following is a brief introduction to the drawings required for the description of the embodiments. Obviously, the drawings described below are only some examples or embodiments of the present application. For ordinary technicians in this field, the present application can also be applied to other similar scenarios based on these drawings without creative work. Unless it is obvious from the language environment or otherwise explained, the same reference numerals in the figures represent the same structure or operation.

[0039] As shown in this application, unless the context clearly indicates an exception, the words "a", "an", "a kind" and / or "the" do not refer to the singular, but also include the plural. Generally speaking, the terms "include" and "comprise" only indicate the inclusion of the steps and elements that have been clearly identified, and these steps and elements do not constitute an exclusive list, and the method or device may also include other steps or elements.

[0040] Unless otherwise specifically stated, the relative arrangement, numerical expressions and numerical values ​​of the parts and steps set forth in these embodiments do not limit the scope of the application. Meanwhile, it should be understood that, for ease of description, the sizes of the various parts shown in the accompanying drawings are not drawn according to actual proportional relationships. The technology, methods and equipment known to those of ordinary skill in the relevant art may not be discussed in detail, but in appropriate cases, the technology, methods and equipment should be considered as a part of the specification. In all examples shown and discussed here, any specific value should be interpreted as being merely exemplary, rather than as a limitation. Therefore, other examples of exemplary embodiments may have different values. It should be noted that similar reference numerals and letters represent similar items in the following drawings, and therefore, once a certain item is defined in an accompanying drawing, it does not need to be further discussed in subsequent drawings.

[0041] In addition, it should be noted that the use of words such as "first" and "second" to define components is only for the convenience of distinguishing the corresponding components. If not otherwise stated, the above words have no special meaning and cannot be understood as limiting the scope of protection of this application. In addition, although the terms used in this application are selected from well-known and commonly used terms, some terms mentioned in the specification of this application may be selected by the applicant at his or her discretion, and their detailed meanings are explained in the relevant parts of the description of this article. In addition, it is required to understand this application not only by the actual terms used, but also by the meaning implied by each term.

[0042] The simulation system for evaluating cardiac ablation schemes of the present application can be used for preoperative evaluation of ablation surgery and guide the formulation of ablation schemes, and can also evaluate the actual rhythm effect of the heart after the ablation surgery has been performed. The simulation system of the present application can be used for simulating arrhythmia diseases such as ventricular and atrial fibrillation and tachycardia, but is not limited thereto.

[0043] Figure 1 is a block diagram of a simulation system for evaluating cardiac ablation schemes in the present application. Figure 1 As shown, the simulation system 100 includes an image acquisition module 110 , a three-dimensional mesh reconstruction module 120 , a fiber rotation direction calculation module 130 , an electrophysiological parameter simulation module 140 and a dynamic simulation module 150 . Among them, the image acquisition module 110 is used to acquire multiple chest images, which include cardiac angiography images; the three-dimensional grid reconstruction module 120 is used to mark the ablation points of the multiple chest images, and obtain the three-dimensional grid of the heart by the finite element method, the ablation points include the estimated ablation points to be ablated, and / or the actual ablation points that have been ablated; the fiber rotation calculation module 130 is used to calculate the myocardial fiber rotation according to the three-dimensional grid of the heart, and obtain the direction vector of the myocardial fiber rotation corresponding to each grid point on the three-dimensional grid of the heart; the electrophysiological parameter simulation module 140 is used to simulate the cardiac action potential excitation process according to the electrophysiological parameters, the geometric position of the grid points, and the fiber rotation corresponding to the grid points, wherein the electrophysiological parameters include cell-level parameters and tissue-level parameters, the cell-level parameters include the current parameters of the ion channel, and the tissue-level parameters include the conductivity of different heart regions, pacing stimulation parameters and the propagation rate of the electrical signal set according to the direction vector of the myocardial fiber rotation; the dynamic simulation module 150 is used to simulate the cardiac electrical signal according to the cardiac action potential excitation process to obtain a dynamic simulation result.

[0044] The following is a detailed description of each module in the simulation system 100.

[0045] The chest images that can be acquired by the image acquisition module 110 refer to images including cardiac angiography images. Preferably, the normal area and the abnormal area can be clearly distinguished from the chest image. The normal area refers to the area including normal cells. The abnormal area is, for example, the area including abnormal conditions such as cell infarction, scar tissue, and fibrotic damage. Multiple chest images means greater than or equal to 2 images. It can be understood that in order to obtain a three-dimensional grid of the heart and establish a three-dimensional model, the number of images should be as large as possible. In order to achieve high-precision simulation results, the original image obtained should have high resolution and high definition.

[0046] This application does not limit how to obtain these chest images.

[0047] In some embodiments, the image acquisition module 110 is specifically implemented as an MRI device, and more specifically, can be implemented as an LGE-MRI device for obtaining LGE-MRI enhanced magnetic resonance images.

[0048] Late Gadolinium Enhancement (LGE)-MRI uses gadolinium-bound macromolecules as contrast agents. Since different areas of the heart accumulate different amounts of contrast agents after using them, the infarcted area can be more clearly distinguished from normal cells. After 10 to 30 minutes of contrast agent use in patients, the contrast agent, relying on the shortening of T1 relaxation time and the different regional distribution patterns of gadolinium-based contrast agents in the extracellular space of myocardial cells, offsets its inherent signal by applying an inversion pulse to the myocardium after a certain period of time, while making the signal of the infarcted tissue brighter. For myocardial infarction, the deposition of collagen increases the extracellular space due to the occurrence of myocardial infarction, so there is a significant difference in the content of gadolinium contrast agents. After weighted imaging, MRI has areas of delayed gadolinium-enhanced signals, revealing areas of damage such as cell infarction, scar tissue, and fibrosis. LGE-MRI can evaluate myocardial scarring and regional myocardial fibrosis. Clinically, LGE-MRI is used to detect and characterize a variety of myocardial diseases, including myocardial infarction, myocarditis, cardiac sarcoidosis, and cardiac amyloidosis.

[0049] The three-dimensional grid reconstruction module 120 obtains the required multiple chest images from the image acquisition module 110. In some embodiments, the multiple chest axial enhanced MRI slice images are segmented into the left ventricular endocardium, the right ventricular endocardium and the epicardium, which is equivalent to marking the boundaries of the left ventricular endocardium, the right ventricular endocardium and the epicardium on each chest image. For other embodiments, for example, for the atrial region, the same operation can also be used to mark the left atrium, the right atrium and the outer surface of the atrium. At the same time, the abnormal areas on the MRI slice images can also be marked by finite element software. The abnormal areas include ablation points, and the ablation points include estimated ablation points to be ablated, and / or actual ablation points where ablation has been performed.

[0050] For preoperative evaluation, the ablation point is the estimated ablation point where ablation will be performed.

[0051] For postoperative evaluation, the ablation points are the actual ablation points where ablation has been performed.

[0052] In some cases, assuming that the electrophysiological signals of a patient who has undergone ablation surgery are abnormal and a second ablation surgery may be required, the ablation points can include both the estimated ablation points for the second ablation surgery and the actual ablation points that have been ablated. In this way, the actual ablation points that have been ablated and the estimated ablation points that will be ablated can be considered simultaneously during simulation, which is closer to the actual situation.

[0053] The present application does not limit how the three-dimensional mesh reconstruction module 120 marks the ablation points. The three-dimensional mesh reconstruction module 120 may provide an algorithm for marking, or provide a graphical user interface for manual marking.

[0054] It should be noted that the three-dimensional mesh reconstruction module 120 obtains a three-dimensional mesh of the heart by a finite element method, which corresponds to the entire heart and includes a plurality of mesh points. The present application does not limit the mesh form, which may be a triangular mesh, a tetrahedral mesh, a hexahedral mesh, etc. The finite element method used is a common method in the art and will not be described in detail here.

[0055] Myocardial fiber handedness is closely related to the propagation of ECG signals. The fiber handedness calculation module 130 is used to calculate the myocardial fiber handedness corresponding to each grid point based on the reconstructed three-dimensional heart grid, and use the direction vector of the myocardial fiber handedness as a parameter of the grid point.

[0056] In one embodiment, the fiber rotation direction calculation module 130 calculates the myocardial fiber rotation direction based on the three-dimensional grid of the heart based on the LD-RB algorithm. Specifically, the rule-based Laplace-Dirichlet algorithm, namely the LD-RB (Laplace-DirichletRuled-based) algorithm, is based on the human myocardial fiber rotation direction data obtained from the cardiac physiological basis research, and obtains the statistical myocardial fiber distribution rules. The parameters of the distribution rules are determined by the three-dimensional grid of the heart, thereby obtaining the fiber rotation direction that meets the relevant rules. This calculation process does not require obtaining special physiological measurement data of the patient.

[0057] The LD-RB algorithm used in this application follows the following fiber handedness distribution rules:

[0058] 1. The longitudinal fibers in the ventricular wall are oriented parallel to the endocardial and epicardial surfaces.

[0059] 2. α is the spiral angle relative to the counterclockwise circumferential direction of the heart when viewed from the base toward the apex. The direction of the longitudinal fibers through the ventricular wall will have a continuous change from the endocardium (-α) to the epicardium (+α). In one embodiment, α is 60°.

[0060] 3. The transverse fiber direction is perpendicular to the longitudinal fiber direction, which is characterized by the angle β, where β is defined as the angle of the ventricle relative to the outer transmural axis.

[0061] 4. The transverse and longitudinal directions of the fiber rotation are perpendicular to the normal of the sheet.

[0062] 5. The fiber rotation direction of the interventricular septum is continuous with the ventricular wall.

[0063] Under the premise of the above five rules, the LDRB algorithm defined above requires four functions as input for the rotation of ventricular muscle fibers, which represent the required α and β angles in the ventricular wall (w) and septum (s), respectively. Here α and β are in degrees, h represents the transmural depth, which is between 0 and 1, with 0 representing the ventricular outer surface and 1 representing the ventricular outer surface. Epi represents the epicardium and Endo represents the endocardium. The four functions are as follows:

[0064] α s (h) = α endo (1-h)-α endo *h (1)

[0065] α w (h) = α endo (1-h)+α epi *h (2)

[0066] β s (h) = β endo (1-h)-β endo *h (3)

[0067] β w (h) = β endo (1-h)+β epi *h (4)

[0068] As boundary conditions, each angle needs to set the value α on the inner and outer membranes w (0) = α s (0) = -α s (1) modulus 180°, and β w (0) = β s (0) = -β s (1) modulo 180°. In the above four functions (1) to (4), except for α and β, the other parameters are known, for example, they can be obtained from the three-dimensional mesh model of the heart. Therefore, α and β can be solved as the rotation angle of the myocardial fiber handedness corresponding to each three-dimensional mesh point.

[0069] The three-dimensional grid of the heart is in the global coordinate system, in which the epicardium, left and right ventricular inner wall, ventricular base and apex are marked. Through coordinate transformation, the direction vector of the myocardial fiber rotation at each grid point can be obtained according to the local coordinate system.

[0070] The electrophysiological parameter simulation module 140 can obtain the required electrophysiological parameters, the geometric position of the grid points, and the fiber rotation direction corresponding to the grid points. For example, the electrophysiological parameter simulation module 140 and the three-dimensional grid reconstruction module 120 and the fiber rotation direction calculation module 130 are all connected in communication, and the geometric position of the grid points can be directly obtained from the three-dimensional grid reconstruction module 120, and the fiber rotation direction corresponding to the grid points can be obtained from the fiber rotation direction calculation module 130. The electrophysiological parameters can be built into the electrophysiological parameter simulation module 140, or set by the user through a graphical user interface.

[0071] In the cellular dimension, electrophysiological parameters include current parameters of ion channels. Specifically, in one embodiment, the electrophysiological parameter simulation module 140 obtains the current parameters of the ion channels corresponding to each grid point on the three-dimensional grid of the heart based on the ion channel model. More specifically, in one embodiment, the ion channel model is the TP06 model. The transmembrane domain of the ion channel constitutes an aqueous pore that allows ions to pass smoothly. It can only allow specific ions to pass through, and controls the transmembrane flow of specific ions by opening and closing a "gate". Many different types of ion channels can usually be found on the myocardial cell membrane, among which sodium (Na+), potassium (K+) and calcium (Ca2+) channels are the main types of ion channels. This application does not limit the type of ion channels. The ion current satisfies the basic Ohm's law, I x =G x p(V m -E x ), where I x represents the ion current of ion x, G x represents the maximum conductivity of ion x, V m represents the transmembrane potential, E x represents the equilibrium potential of ion x, and p represents the open probability coefficient of the ion channel. Ion channels are voltage-controlled, that is, the opening and closing of ion channels are affected by the transmembrane potential on both sides of the membrane.

[0072] In the tissue dimension, the electrophysiological parameters include the conductivity of different cardiac regions. It should be noted that the conductivity of myocardial cells in normal areas is different from that in abnormal areas. In one embodiment, the conductivity of different areas in the three-dimensional grid of the heart can be determined based on the segmented chest image, wherein the conductivity corresponding to the ablation point is different from the conductivity corresponding to the normal myocardial cells. As mentioned above, the infarct area and normal cells can be clearly distinguished in the LGE-MRI image. The infarct area can be identified through the images in turn. The fibrotic scar in the infarct area is mainly composed of an extracellular matrix composed of collagen fibers, which is not active in electrophysiological properties, so it is usually regarded as an insulator, that is, the conductivity is approximately 0. In order to simplify the calculation process, the conductivity at the ablation point is set to 0 in this application. Correspondingly, the conductivity of normal cells can also be set to an equal value, such as the average value of the actual conductivity.

[0073] In the edge area of ​​fibrotic matrix scar, cells are still active, but the action potential of the cardiomyocytes in the edge area will be different from that of normal cardiomyocytes. The speed and height of the potential rise during depolarization will be reduced, and the APD (Action Potential Duration) time will be prolonged. The regional points caused by ablation will also be regarded as insulating areas, so the setting of electrophysiological signal parameters is also applicable to the scar itself or the ablation results.

[0074] In the tissue dimension, the electrophysiological parameters also include the propagation rate of the electrical signal set according to the direction vector of the myocardial fiber handedness calculated by the fiber handedness calculation module 130. In some embodiments, the propagation rate is also set with an average value. For different grid points, due to the anisotropic characteristics of the fiber handedness, the propagation rates corresponding to different grid points are the same but in different directions.

[0075] In the organizational dimension, electrophysiological parameters also include pacing stimulation parameters. Specifically, pacing stimulation parameters can be set according to the clinical requirements of the patient. Pacing stimulation parameters vary according to the stimulation scheme. For programmed action potential stimulation schemes, you can consider trying multiple designs and perform twin replication of the patient's condition in the simulation. In order to achieve the purpose of personalization for the patient, several different types of designs may be used:

[0076] 1. Graded incremental stimulation: pacing at a frequency 10 to 20 beats / min faster than the subject's basic heart rate, pacing for 30 to 60 seconds each time. After an interval of 1 to 2 minutes, the frequency increases by 10 to 20 beats / min to enter the next cycle until a certain part of the heart reaches the effective refractory period or the longest sinus node recovery time occurs.

[0077] 2. Continuous incremental stimulation: Use a lower frequency to pace, then gradually increase the pacing frequency, maintain 1:1 pacing capture, until a certain part of the heart reaches an effective refractory period. This stimulation method can reduce the impact of overspeed inhibition and is often used for refractory period determination.

[0078] 3. Short burst rapid stimulation: Rapid pacing 10 to 20 times at a frequency 30 beats / minute faster than the heart rate of the person being examined to terminate rapid arrhythmia.

[0079] 4. Train stimulation: consists of several groups of short, rapid stimulation pulses.

[0080] After considering the above different types of design factors, based on the doctor's consensus, a programmed premature stimulation scheme can be used to try to induce the corresponding condition during the simulation process, in order to achieve the purpose of patient personalization. Programmed premature stimulation is to isolate single or multiple premature (premature) stimulations in the native heart rhythm or the basic paced heart rhythm.

[0081] The pre-procedural stimulation scheme adopted in this application is configured to include the following steps:

[0082] Step 1: Select the location of the stimulation point;

[0083] Step 2: Use S1S2 stimulation. The stimulation sequence includes 6 consecutive stimulations of the stimulation point by a pacing stimulation S1 with a cycle length of T1. After the last S1 stimulation, apply an S2 stimulation. The cycle length of the S2 stimulation is T2. If no abnormality is induced after the S2 stimulation, the cycle length T2 of the S2 stimulation is decreased according to a preset step length, and the S2 stimulation is applied cyclically until an abnormality is induced. In some cases, T1=600ms, the initial value of T2 is 250ms, and the decreasing step length is 10ms.

[0084] Step 3: If the S2 stimulation still fails to induce abnormality, then apply the S3 stimulation, and decrease the cycle length T3 of the S2 stimulation according to the preset step length, and cyclically apply the S3 stimulation until the abnormality is induced.

[0085] The reason why abnormality cannot be induced may be that the myocardial tissue at the stimulation point is still in the refractory period. In some cases, a threshold is set for T2 of S2. When abnormality cannot be induced even after reaching the threshold, S2 stimulation is stopped and S3 stimulation is applied instead. For example, when T2 decreases to <150ms, S2 stimulation is stopped.

[0086] Step 4: If the S3 stimulation still fails to induce abnormality, then the S4 stimulation is applied, and the cycle length T4 of the S4 stimulation is decreased according to the preset step length, and the S4 stimulation is applied cyclically until the abnormality is induced;

[0087] Step 5: If S4 stimulation still fails to induce abnormality, it is considered that the stimulation point cannot induce the desired arrhythmia. In this case, the stimulation point can be changed and a new round of stimulation can be performed.

[0088] This application does not impose any restrictions on the location of the stimulation point, the intensity of stimulation, etc., which can be determined by the doctor based on the standards.

[0089] Based on the above content, the electrophysiological parameter simulation module 140 can simulate the cardiac action potential excitation process and program the cardiac action potential excitation process. On this basis, the dynamic simulation module 150 can simulate the cardiac electrical signal according to the cardiac action potential excitation process to obtain the simulation result. The dynamic simulation module 150 can visualize the simulation result.

[0090] Specifically, in one embodiment, the dynamic simulation module 150 simulates the cardiac electrical signal based on the Reaction-Diffusion diffusion process, and the simulation process includes:

[0091] Step 10: The transmembrane position and reset voltage on the time-space evolution scale are used as contraction variables u, v and nonlinear response functions f, g to satisfy equations (5) and (6):

[0092]

[0093] Among them, D1 and D2 are the coefficients corresponding to the shrinkage variables u and v, is the Laplace operator, More specifically, D1 and D2 are used to characterize the conductivity of the corresponding diffusion process in the transmembrane current and slow local recovery current in the ion channel model.

[0094] It should be noted that cardiac electrical activity in a single cell involves the coordinated transport of a large number of ions through various biological channels. Electrical waves also propagate from cells to the entire heart, combining excitation with contraction and blood expulsion. Some simplifying assumptions and rules, such as cellular automaton models, limit their ability to simulate cardiac electrical activity, especially at the level of the entire heart. The reaction-diffusion model describes the distribution changes of dynamic variables in space and time, which includes two processes: (1) reaction process: dynamic variables interact to transform, and (2) diffusion process: dynamic variables diffuse in space. The potential propagation process determined by the reaction-diffusion model is reflected in the time-space evolution as a continuously diffusing contour map boundary. Therefore, the reaction-diffusion model can simulate the electrical activity of the whole heart more realistically. However, the reaction-diffusion model is more difficult to implement on complex geometric shapes and has higher computational costs. The anisotropy introduced by the fiber handedness in this application does not increase the complexity of the reaction-diffusion model and can make the calculation results more realistic.

[0095] Step 20: Discretize equations (5) and (6) in step 10 in the time dimension, set t = 1, 2, ..., T, and record the grid point numbers i = 1, 2, ..., N in the three-dimensional grid of the heart, and obtain the following equations:

[0096]

[0097] Step 30: Write the equation in step 20 as the following matrix:

[0098]

[0099] in, Similarly, the matrix forms of B1, B2, F, and G can be determined. B1 and B2 are the corresponding coefficient matrices. The specific parameters of this matrix need to be determined by the finite element form of the Laplacian operator on the corresponding mesh. For a triangular mesh, the Laplacian operator acts on a variable u. where Δ is M is called the Mass matrix, and L is the Laplacian matrix.

[0100] where θ ij and The mass matrix M is a diagonal matrix related to the area of ​​the finite element triangle where each grid point is located. in To find the area of ​​the triangle, we use the cross product of the measured coordinates. After calculating the fiber spin, in order to add the anisotropy of ECG signal propagation, we can use the fiber spin as the local coordinate basis vector. In this way, the local coordinates will need to be obtained by multiplying the original measured coordinates by the metric. The basis vector It is longitudinal along the fiber direction. is the transverse direction of the fiber, It is the perpendicular normal to the plane formed by the direction vectors of the two fiber rotation directions. Therefore, the vector represented by the calibration coordinates will be expressed in local coordinates as At this time, the coefficients of each item in formula (9) can be determined and satisfy the form of Ax=b, so the solution can be obtained.

[0101] Step 40: Use optimization method to solve the matrix in step 30 and obtain U t and V t .

[0102] The solution can be obtained by using an optimization method commonly used in the art, and this application does not limit the specific optimization method.

[0103] In one embodiment, in order to simulate the dynamic process in time and space, a nonlinear reaction-diffusion model is used. Specifically, in step 10, the Fitzhugh-Nagumo model is used to change the functions f and g (giving the functions f and g a specific form) to obtain:

[0104] f(u,v)=C1u(1-u)(ua)-C2uv (10)

[0105] g(u,v)=b(u-dv) (11)

[0106] Formula (10) and formula (11) correspond to functions f and g, respectively. Among them, C1 and C2 are parameters related to the curve shape of function f, 1 / b is the average time constant of the slow current passing through the ion channel, a is a parameter related to the resting voltage, and d is the ratio of the transmembrane current and the slow local recovery current in the ion channel. Parameters a, b, C1, C2, and d are all related to the specific ion channel model, and their actual values ​​need to be determined by the shape of the action potential curve of the model. According to the functions f and g obtained from the above formulas (10) and (11), the diffusion equations of the variables u and v shown in equations (5) and (6) satisfy:

[0107]

[0108] In an embodiment using the TP06 ion channel model, the correlation coefficients may be set as follows: D1=1, D2=0, a=0.13, b=0.013, C1=0.26, C2=0.1, d=1.

[0109] The fiber handedness will determine the Laplace operator The numerical solution at each step of discretization, and then at each step of simulation, reflects the anisotropy of the electrical signal propagating along the fiber rotation. It should be noted that in ordinary diffusion equations, anisotropy is brought about by conductivity D1 as a tensor matrix. This application resorts to the Laplace operator instead of the conductivity tensor for anisotropy. The fiber rotation will determine the direction of the nearest neighbor point of each spatial grid point in the spatial discretization, and change the metric of the local coordinates to change the anisotropy of the Laplace operator in the numerical calculation. This method of introducing anisotropy does not increase the computational difficulty of the reaction-diffusion equations (5) and (6). Furthermore, considering the anisotropy of the electrical signal speed in different directions, the anisotropy brought about by the fiber rotation also reflects the propagation speed of the electrical signal between each pair of nearest neighbor points. Therefore, the simulation results obtained by the simulation system of this application are more likely to be close to the propagation process of the actual electrophysiological signal, improve the authenticity of the simulation results, and can be more effectively applied to clinical electrophysiological judgment.

[0110] like Figure 1 As shown, in some embodiments, the simulation system 100 of the present application may further include an evaluation module 160. The evaluation module 160 is used to evaluate the cardiac ablation scheme based on the dynamic simulation results, and the evaluation results include whether to change the ablation point. Specifically, for patients who plan to undergo ablation, the simulation system 100 can obtain the patient's chest image and obtain the dynamic simulation results for the patient. The evaluation module 160 can compare the dynamic simulation results with the normal cardiac electrophysiological process, and evaluate the actual impact of the lesion point or the estimated ablation point of the ablation scheme, thereby assisting the doctor in determining whether the ablation point needs to be changed. If a change is required, the evaluation module 160 transmits the evaluation results to the three-dimensional mesh reconstruction module 120, and the three-dimensional mesh reconstruction module 120 needs to re-mark the ablation point, and then the fiber rotation direction calculation module 130, the electrophysiological parameter simulation module 140 and the dynamic simulation module 150 simulate the changed ablation point in turn to obtain a new simulation result, which is then evaluated.

[0111] Figure 2 FIG. 1 is an exemplary flow chart of using a simulation system according to an embodiment of the present application for evaluating a cardiac ablation plan. Figure 2 As shown, the process includes the following steps:

[0112] Step S210: Obtaining estimated ablation points.

[0113] Step S220: obtaining actual ablation points.

[0114] It should be noted that step S210 and step S220 are performed according to actual needs. If the patient to be evaluated is going to undergo an ablation surgery, the estimated ablation points in the ablation plan are obtained and the ablation plan is evaluated. If the patient to be evaluated has already undergone an ablation surgery, the actual ablation points are obtained and the ablation plan that has been implemented is evaluated. For patients who need a second ablation, the estimated ablation points and the actual ablation points can be obtained at the same time.

[0115] Step S230: using the simulation system 100 described above to perform simulation.

[0116] Step S240: Evaluate according to the simulation results. For the evaluation of the estimated solution, if the evaluation result indicates that the ablation point needs to be modified, the estimated solution is modified before simulation.

[0117] Step S250: If the simulation result indicates that the ablation plan is feasible, a personalized ablation plan for the patient can be obtained.

[0118] Figure 2 The figure is only one application mode of the simulation system 100 of the present application, and is not used to limit the specific application scenario of the simulation system 100.

[0119] The simulation system 100 of the present application constructs a three-dimensional grid of the heart through images, calculates the rotation direction of myocardial fibers, and obtains a cardiac action potential excitation process with fiber rotation anisotropy, thereby introducing anisotropic cardiac potential diffusion processes in the RD diffusion model for simulation calculations, simulating the actual electrophysiological properties of myocardial fibers, and greatly improving the numerical accuracy. At the same time, it also avoids the instability of the results obtained by the boundary element method, circumvents the defects of the ECG inverse problem, and can simulate the characteristics of actual cases more finely. In addition, in the RD diffusion model, known physiological parameters can be conveniently applied, and relevant parameters can be flexibly adjusted in the simulation calculation, which can more realistically simulate the evolution of the target patient's disease. The present application applies the anisotropy of fiber rotation to the relevant parameters of the diffusion model, which improves the authenticity of the simulation, so that the simulation results are not only qualitatively displayed, but also have the accuracy of digital twin quantitative research on diseases.

[0120] The basic concepts have been described above. Obviously, for those skilled in the art, the above application disclosure is only an example and does not constitute a limitation of the present application. Although not explicitly stated herein, those skilled in the art may make various modifications, improvements and amendments to the present application. Such modifications, improvements and amendments are suggested in the present application, so such modifications, improvements and amendments still belong to the spirit and scope of the exemplary embodiments of the present application.

[0121] At the same time, the present application uses specific words to describe the embodiments of the present application. For example, "one embodiment", "an embodiment", and / or "some embodiments" refer to a certain feature, structure or characteristic related to at least one embodiment of the present application. Therefore, it should be emphasized and noted that "one embodiment" or "an embodiment" or "an alternative embodiment" mentioned twice or more in different positions in this specification does not necessarily refer to the same embodiment. In addition, some features, structures or characteristics in one or more embodiments of the present application can be appropriately combined.

[0122] Some aspects of the present application may be performed entirely by hardware, entirely by software (including firmware, resident software, microcode, etc.), or by a combination of hardware and software. The above hardware or software may be referred to as "data blocks", "modules", "engines", "units", "components" or "systems". The processor may be one or more application specific integrated circuits (ASICs), digital signal processors (DSPs), digital signal processing devices (DAPDs), programmable logic devices (PLDs), field programmable gate arrays (FPGAs), processors, controllers, microcontrollers, microprocessors, or combinations thereof. In addition, various aspects of the present application may be expressed as computer products located in one or more computer-readable media, which include computer-readable program codes. For example, computer-readable media may include, but are not limited to, magnetic storage devices (e.g., hard disks, floppy disks, tapes ...), optical disks (e.g., compact disks CDs, digital versatile disks DVDs ...), smart cards, and flash memory devices (e.g., cards, sticks, key drives ...).

[0123] A computer-readable medium may include a propagated data signal containing computer program code, such as in baseband or as part of a carrier wave. The propagated signal may have a variety of manifestations, including electromagnetic, optical, etc., or a suitable combination. A computer-readable medium may be any computer-readable medium other than a computer-readable storage medium, which may be connected to an instruction execution system, device or apparatus to communicate, propagate or transmit a program for use. The program code on the computer-readable medium may be propagated via any suitable medium, including radio, cable, fiber optic cable, radio frequency signal, or similar medium, or any combination of the above mediums.

[0124] Similarly, it should be noted that in order to simplify the description of the disclosure of this application and thus help understand one or more embodiments of the invention, in the above description of the embodiments of the present application, multiple features are sometimes combined into one embodiment, figure or description thereof. However, this disclosure method does not mean that the subject of the present application requires more features than the features mentioned. In fact, the features of the embodiment are less than all the features of the single embodiment disclosed above.

[0125] In some embodiments, numbers describing the number of components and attributes are used. It should be understood that such numbers used in the description of the embodiments are modified by the modifiers "about", "approximately" or "substantially" in some examples. Unless otherwise specified, "about", "approximately" or "substantially" indicate that the numbers are allowed to vary by ±20%. Accordingly, in some embodiments, the numerical parameters used in this application are approximate values, which may change according to the required characteristics of individual embodiments. In some embodiments, the numerical parameters should take into account the specified significant digits and adopt the general method of retaining the digits. Although the numerical domains and parameters used to confirm the breadth of the range in some embodiments of the present application are approximate values, in specific embodiments, the setting of such numerical values ​​is as accurate as possible within the feasible range.

[0126] Although the present application has been described with reference to the current specific embodiments, ordinary technicians in this technical field should recognize that the above embodiments are only used to illustrate the present application, and various equivalent changes or substitutions may be made without departing from the spirit of the present application. Therefore, as long as the changes and modifications to the above embodiments are within the essential spirit of the present application, they will fall within the scope of the present application.

Claims

1. A simulation system for evaluating cardiac ablation schemes, characterized in that: include: An image acquisition module, used for acquiring a plurality of chest images, wherein the chest images include cardiac angiography images; A three-dimensional mesh reconstruction module, used for marking ablation points on the multiple chest images and obtaining a three-dimensional mesh of the heart by a finite element method, wherein the ablation points include estimated ablation points to be ablated and / or actual ablation points where ablation has been performed; A fiber handedness calculation module, used to calculate the handedness of myocardial fibers according to the three-dimensional cardiac grid, and obtain a direction vector of the handedness of myocardial fibers corresponding to each grid point on the three-dimensional cardiac grid; An electrophysiological parameter simulation module, used to simulate the cardiac action potential excitation process according to the electrophysiological parameters, the geometric positions of the grid points, and the fiber rotation directions corresponding to the grid points, wherein the electrophysiological parameters include cell-level parameters and tissue-level parameters, the cell-level parameters include the current parameters of the ion channels, and the tissue-level parameters include the conductivity of different cardiac regions, pacing stimulation parameters, and the propagation rate of the electrical signal set according to the direction vector of the myocardial fiber rotation direction; The dynamic simulation module is used to simulate the cardiac electrical signal according to the cardiac action potential excitation process to obtain a dynamic simulation result.

2. The simulation system according to claim 1, characterized in that: The chest image includes an LGE-MRI enhanced image.

3. The simulation system according to claim 1, characterized in that: The three-dimensional grid reconstruction module is also used to segment the multiple chest images, and the segmented boundaries include the left ventricular endocardium, the right ventricular endocardium and the epicardium.

4. The simulation system according to claim 1, characterized in that: The fiber handedness calculation module calculates the myocardial fiber handedness according to the three-dimensional heart grid based on the LD-RB algorithm.

5. The simulation system according to claim 1, characterized in that: The electrophysiological parameter simulation module obtains the current parameter of the ion channel corresponding to each grid point on the three-dimensional grid of the heart based on the ion channel model.

6. The simulation system according to claim 3, characterized in that: The electrophysiological parameter simulation module determines the conductivity of different regions in the three-dimensional grid of the heart according to the segmented chest image, wherein the conductivity corresponding to the ablation point is different from the conductivity corresponding to normal myocardial cells.

7. The simulation system according to claim 6, characterized in that: The electrical conductivity corresponding to the ablation point is 0.

8. The simulation system according to claim 1, characterized in that: The pacing stimulation parameters include parameters corresponding to a pre-programmed stimulation scheme, and the pre-programmed stimulation scheme is configured to include: Step 1: Select the location of the stimulation point; Step 2: Use S1S2 stimulation, the stimulation sequence includes 6 consecutive stimulations of the stimulation point by the pacing stimulation S1 with a cycle length of T1, and after the last S1 stimulation, apply an S2 stimulation, the cycle length of the S2 stimulation is T2, if no abnormality is induced after the S2 stimulation, then the cycle length T2 of the S2 stimulation is decreased according to a preset step length, and the S2 stimulation is applied cyclically until the abnormality is induced; Step 3: If the S2 stimulation still fails to induce abnormality, then the S3 stimulation is applied, and the cycle length T3 of the S2 stimulation is decreased according to the preset step length, and the S3 stimulation is applied cyclically until the abnormality is induced; Step 4: If the S3 stimulation still fails to induce abnormality, then the S4 stimulation is applied, and the cycle length T4 of the S4 stimulation is decreased according to the preset step length, and the S4 stimulation is applied cyclically until the abnormality is induced; Step 5: If S4 stimulation still fails to induce abnormalities, it is considered that the stimulation point cannot induce the desired arrhythmia.

9. The simulation system according to claim 1, characterized in that: The dynamic simulation module simulates the cardiac electrical signal based on the Reaction-Diffusion diffusion process, and the simulation process includes: Step 10: The transmembrane position and reset voltage on the time-space evolution scale are used as contraction variables u, v and nonlinear response functions f, g to satisfy the equation: Among them, D1 and D2 are the coefficients corresponding to the shrinkage variables u and v, is the Laplace operator, Step 20: Discretize the equation in step 10 in the time dimension, set t = 1, 2, ..., T, and record the grid point sequence i = 1, 2, ..., N in the three-dimensional grid of the heart, and obtain the following equation: Step 30: Write the equation in step 20 as the following matrix: in, Step 40: Use optimization method to solve the matrix in step 30 and obtain U t and V t .

10. The simulation system according to claim 9, characterized in that: In step 10, the Fitzhugh-Nagumo model is adopted so that: f(u,v)=C1u(1-u)(ua)-C2uv g(u,v)=b(u-dv) Among them, C1 and C2 are parameters related to the curve shape of the function f, 1 / b is the average time constant of the slow current passing through the ion channel, a is a parameter related to the resting voltage, and d is the ratio of the transmembrane current to the slow local recovery current in the ion channel.

11. The simulation system according to claim 1, characterized in that: Also includes: An evaluation module is used to evaluate the cardiac ablation plan according to the dynamic simulation result, and the evaluation result includes whether to change the ablation point.

Images