Patent Application
20250169782
The invention pertains to myocardial perfusion imaging (MPI) and poromechnics. In particular, the invention involves the use of non-invasive measurements (for example, computed angiography or a nuclear magnetic resonance) to build a model comprising a porous domain of a human heart to be used in computational fluid dynamics. The results of the invention allow to diagnose the human heart. For example, to diagnose a coronary heart disease.
One of the major challenges for the application of perfusion assessment methods in medical diagnosis is a coronary heart disease (CHD). According to the National Center for Health Statistics, 5,6% of U.S. adults had a coronary heart disease in 2018. Prevalence increased with age and almost 24% of people at the age of 75 and over had a CHD. The disease affected men more often than women—7.4% vs 4.1% of all U.S. adults (according to Villarroel M A, Blackwell D L, Jen A (2019) Tables of Summary Health Statistics for U.S. Adults: 2018 National Health Interview Survey, National Center for Health Statistics. Accessed 15 Dec. 2021). In 2018, a CHD caused 365 744 deaths in the U.S. which accounts for almost 13% of all deaths therein (see Centers for Disease Control and Prevention, National Center for Health Statistics, National Vital Statistics System: public use data file documentation: mortality multiple cause-of-death micro-data files. Accessed 15 Dec. 2021). In 2017, 623 276 people died from an ischemic heart disease in the European Union which accounts for almost 12% of all deaths therein in 2018 (see Eurostat (2021) Causes of death-deaths by country of residence and occurrence [HLTH_CD_ARO]. Accessed 15 Dec. 2021).
A heart is the center of the circulatory system. A human heart beats about 100 000 times per day which sums to about 3 billion beats in a lifetime. Pumping of a heart and a heartbeat are caused by alternating contractions and relaxations of the myocardium. This pumping ensures continuity of blood movement in the vessels. Thus, a heart is responsible for fulfilling various functions in a human body: provision of oxygen and nutrients, removal of metabolites, protection from foreign bodies, toxic substances and against exsanguination. Moreover, a heart is responsible for maintenance of constant physical and chemical conditions in a human body, especially temperature and pH (according to Barrett K E et all (2019) Ganong's Review of Medical Physiology, 26th Edition, McGraw-Hill Education). Indeed, the modern medical diagnostics offer a wide spectrum of methods and parameters for a functional assessment of a heart. Despite this, like a hundred and more years ago, the key role in medical diagnostics is played by a parameter that describes the amount of blood pumped within a certain time interval, usually within a single cycle or a minute. As long as this parameter is maintained at an adequate level, the prerequisite for other organs to maintain their functions is fulfilled. The awareness of this principle and attempts to evaluate this parameter can be found in very early experimental works on physiology. Already in 1761, Albrecht von Haller was studying a pulmonary circulation by an injection of a colored liquid into the vena cava of a fresh animal corpses (see Haller A (1761) Elementa physiologiae corporis humani, Francisci Grasser & Sociorum, Lausanne). The turning point and the foundation of the modern methods of flow diagnostics with indicator-dilution techniques, i.e., the techniques commonly used in modern times to assess perfusion, was the research done by George Stewart (see Argueta E E, Paniagua D (2019) Thermodilution cardiac output, Cardiology in Review, 27(3), p. 138-144) and then, less than a century later, the work done by Kenneth Zierler (see Meier P, Zierler K L (1954) On the theory of the indicator-dilution method for measurement of blood flow and volume, The Journal of Applied Physiology, 6(12), p. 731-744, Zierler K (2000) Indicator dilution methods for measuring blood flow, volume, and other properties of biological systems: a brief history and memoir, Annals of Biomedical Engineering, 28(8), p. 836-848).
A myocardial prefusion, i.e., a blood flow through a specific mass of a tissue, is the basic physiological parameter determined on the basis of cardiological image data. Perfusion imaging—including myocardial perfusion imaging (MPI)—involves injecting a contrast into the blood and recording distribution of changes in its concentration in the cardiovascular system over time. A key issue in a perfusion analysis is the calculation of so-called perfusion parameters or perfusion metrics such as myocardial blood flow MBF, myocardial blood volume MBV, mean transit time MTT and time to peak TTP (following Taylor S H (1966) Measurement of the cardiac output in man, Proc R Soc Med, 59 Suppl(Suppl 1), p. 35-53, Axel L (1980) Cerebral blood flow determination by rapid-sequence computed tomography: theoretical analysis, Radiology, 137(3), p. 679-86, Lee T-Y (2002) Functional CT: Physiological models, Trends Biotechnol, 20, p. S3-S10, Tomandl B F et al. (2003) Comprehensive imaging of ischemic stroke with multisection CT, Radiographics, 23(3), p. 565-92, Konstas A A et al. (2009) Theoretic basis and technical implementations of CT perfusion in acute ischemic stroke. Part 1: Theoretic basis, AJNR Am J Neuroradiol, 30(4), p. 662-8). Let's assume a certain control volume of interest (VOI) with a uniquely defined inlet and an outlet. If a certain amount (m) of a contrast agent was introduce at the initial moment t0=0, and a sample was taken on the outlet side, after the time Δti, in which the concentration c(ti), ti∈Δti of the contrast agent was measured, then the mass of the contrast agent that left VOI is expressed by
where ΔV(ti)=q(ti)·Δti is a volume of blood that flowed through an outlet section. The total mass of the contrast agent that will reach the outlet, after sampling at the outlet and determining the concentration for a sufficiently long time, is expressed by
If samples were taken frequently enough than we have
Hence, based on the mean value theorem we can arrive at
where
that describes the volume (e.g., in liters) of blood being pumped by a heart per unit time (e.g., per minute) as the amount of contrast agent at the inlet divided by the area under concentration curve (AUC) recorded at the outlet.
If we analyze an in vivo recording of the arterial input function (AIF) and the time attenuation curve (TAC): a change of contrast agent quantity that occurs in a certain time interval is directly proportional to the difference between inflow and outflow quantity. This can be expressed by
which determines the balance of the mass flux of a contrast agent in the control volume Vt of the tissue, where the average flow is Qt. The concentration c(t) in the region monitored by the arterial input AIF curve was denoted as ca(t) and the tissue attenuation TAC curve as ct(t), and the concentration in venous outflow as cv(t). Given that relative unit values rather than absolute values are used in the perfusion analysis, it is appropriate to express the mean flow Q accordingly by the mean volume of flowing blood (MBF—myocardial blood flow) in a bulk of a tissue:
where ρt is tissue density. Thus, the output balance equation takes the form of
When analyzing the shape of TAC, it can be noticed that in the initial phase there is a gradual accumulation of the contrast agent in the control volume and the outflow is negligible, and the concentration curve has a strictly increasing character. For this phase, assuming that the analysis is only limited to the part where the outflow is negligibly small, i.e., cv(t)≈0, the balance equation takes the form of
and this assumption can be additionally strengthened by limiting our considerations only to the maximum values
So, if we want to assess perfusion (MBF—myocardial blood flow) in a practical way, we can use the maximum slope method. Then MBF is expressed as the ratio of the maximum slope of TAC for the accumulation phase to the maximum value recorded on AIF, multiplied by the reciprocal of the tissue density ρt
By the same reasoning and by starting from the mass balance and not its flux as in equation (17), we obtain the relation for the unit volume of the vascular bed, which is called the myocardial blood volume (MBV)
Calculation of the perfusion time parameters, e.g., mean MTT and max TTP, is a purely arithmetic operation
TTP is simply a coordinate on the time axis for which the following relation
is satisfied. The maximum slope method is a method with extremely high efficiency and moderate hardware requirements. It is possible to apply it effectively even to very large size data sets. This is of considerable importance for analyzing of computed tomography (CT) perfusion imaging data.
The field of noninvasive clinical cardiac imaging is rapidly evolving. The evolution is focused on quantitative perfusion assessment technologies which can be used for assessing a coronary flow and a myocardial ischaemia. The main methods for myocardial perfusion assessment are the single photon emission computed tomography (SPECT), the positron emission tomography (PET) and the stress dynamic computed tomography perfusion (CTP) (see AlJaroudi W A, Hage F G (2020) Review of cardiovascular imaging in the Journal of Nuclear Cardiology 2019: Positron emission tomography, computed tomography and magnetic resonance, J Nucl Cardiol, 27(3), p. 921-930).
For decades, SPECT has been an established method for clinical cardiovascular imaging. The evaluation of a myocardial ischaemia extent based on SPECT, taken together with measurements of a left ventricular function and volumes, is considered to be an approved method for diagnosis of a myocardial ischaemia. The main goals of the SPECT technical evolution were maximizing sensitivity, reducing of acquisition time and reducing amounts of injected radiopharmaceuticals. A SPECT temporal resolution (which is a frame acquisition time) is equal to 10 s and a voxel size is 10×10×10 mm3. In addition, there was the problem of increasing BMI (body mass index) in patient population. The solution to the problem was cadmium-zinc-telluride solid-state detector technology (CZT) which provided an increased sensitivity, an increased spatial and energy resolution and a decreased acquisition time and with smaller doses of radiopharmaceuticals. CZT SPECT enabled a dynamic imaging and kinetics assessment of perfusion tracers in blood and in the myocardium. With this information, a compartmental modelling and delineation of absolute MBF and a coronary flow reserve are possible (see Dewey M (2020) Clinical quantitative cardiac imaging for the assessment of myocardial ischaemia, Nature Reviews Cardiology, 17(7), p. 427-450).
Another technique for confirming or excluding ischaemia is the PET. PET allows for an accurate in vivo measurement of concentrations of radiopharmaceuticals. A PET temporal resolution (which is a frame acquisition time) is in the range of 1-5 s and a voxel size is 4×4×4 mm3. Different markers are used in examination, but 15O-water is considered to be a clinical gold standard. Modelling of 15O-water kinetics is well-established and tracked concentration thereof is directly proportional to measured image signal. Because of short half-lives of tracers, rest-stress protocols can be completed in 30 min, therefore a total radiation dose is lower than that in SPECT. Images of perfusion defects are sometimes grainy and difficult to interpret. MBF values can be provided by modelling tracer kinetics and correction for limited extraction. Recent advances in PET allowed for computation of parametric images showing voxel-level MBF and facilitated 15O-water PET use in the clinical environment (see Dewey M (2020) Clinical quantitative cardiac imaging for the assessment of myocardial ischaemia, Nature Reviews Cardiology, 17(7), p. 427-450).
The last of these three techniques is CTP. CTP presently has two approaches: static and dynamic. The assessment of myocardial enhancement acquired at a single time point during the first-pass of a contrast sample is referred to as a static CTP imaging. The dynamic CTP, on the other hand, assesses a myocardial perfusion by measuring myocardial enhancement at numerous time intervals during the initial pass of a contrast sample, which is bolus timing resilient. The dynamic CTP allows for a fully quantitative examination of a myocardial perfusion. Due to the recent advancements in a computed tomography (CT) technology, the dynamic CTP of the whole myocardium may now be conducted with a significantly less radiation exposure (see Kitagawa K (2018) Dynamic CT perfusion imaging: state of the art, Cardiovascular Imaging Asia, 2(2), p. 38-48).
All the above-mentioned methods have both advantages and disadvantages. In terms of spatial and temporal resolutions, a quality of an image and an accuracy of diagnosis, PET outperforms SPECT. PET can be also performed in patients with implantable cardioverter defibrillators and in those with pacemakers. However, PET involves high costs, a high radiation exposure and is not so widely available. SPECT is more suitable for a clinical routine, because of radionuclides being less expensive, easier to prepare and having longer half-lives. In addition, SPECT performs better in detection of an ischaemia due to its high sensitivity and specificity. SPECT allows for the assessment of a left ventricular function. The limitations of SPECT are high costs, time consumption, a radiation exposure, limited anatomical information because of a poor spatial resolution, attenuation artefacts producing false positive results (especially in obese patients), possibility of underestimating the actual extent of an ischaemia in patients with a multivessel disease. Unlike other mentioned methods, a CTP imaging is the only non-invasive method for identifying functional significance of coronary stenosis and quantifying it. A single CTP examination provides an integrated anatomic and a functional evaluation. Moreover, it is very fast and a widely available method with a high sensitivity and specificity. CTP allows for detection of smaller perfusion defects due to a high spatial resolution (approximately 0.3 mm) thereof. However, this method suffers from radiation exposure and motion artefacts originating from breathing (breath artefacts), beam hardening and high heart rate artifacts (following Xu C, Yi Y & Wang Y (2020) Clinical Applications of CT Myocardial Perfusion Imaging, Cardiovascular Imaging Asia, 4(4), p. 86-90, Seitun S (2018) CT myocardial perfusion imaging: a new frontier in cardiac imaging, BioMed research international, 2018).
A progress in a computed tomography technology in recent years has allowed for the clinical application of a myocardial perfusion imaging by the computed tomography (CTP). In line with PET and the magnetic resonance imagining (MRI), CTP may provide qualitative or quantitative perfusion data. An assessment of a myocardial perfusion by CT can be achieved through a static or a dynamic scan acquisition. Although both types of acquisition have proven to be clinically applicable, important barriers must be overcome before they will be used routinely. For example, a success of a perfusion assessment using the static CTP is highly dependent on a contrast material bolus timing and motion. In contrast, the dynamic CTP of a myocardial perfusion assessment allows for a fully quantitative analysis of a myocardial perfusion. A myocardium is enhanced at a plurality of time points of first-pass of a contrast agent, in that way a myocardium perfusion is assessed. An assessment of a myocardial perfusion may also be affected by performing stress CTP studies with and without a beta-blockade (see Seitun S (2016) Stress computed tomography myocardial perfusion imaging: a new topic in cardiology, Revista Española de Cardiología (English Edition), 69(2), p. 188-200).
During perfusion examination, difficulties can arise in establishing clear cut-off MBF values. Among researchers, the cut-off values of the MBF derived from CT vary between 75 and 164 mL/100 mL/min. Differences in a study design, a patient coronary risk profile, sample sizes, an FFR threshold value, and prevalence of a CAD may be the reason for this wide range of cut-off MBF values (see Kitagawa K (2018) Dynamic CT perfusion imaging: state of the art, Cardiovascular Imaging Asia, 2(2), p. 38-48). When MBF is uniformly high (i.e., >110 mL/100 mL/min for the dual-source CT and a maximum slope technique), a stress myocardial perfusion can be regarded to be normal. A low MBF (75 mL/100 mL/min) is uncommon, however it is most likely related to a lack of or to an inadequate hyperemic response to vasodilators, which might make ischemia identification difficult. Once sufficient hyperemia has been proven by a high MBF in a remote myocardium, a hemodynamic significance of a coronary artery stenosis can be determined using a relative MBF value in the vascular region. Because of normal MBF patient-to-patient variations, a generally applicable absolute MBF threshold for a myocardial ischemia may not exist (see Kitagawa K (2018) Dynamic CT perfusion imaging: State of the art, Cardiovasc Imaging Asia, 2(2), p. 38-48).
With a CTP imaging, we can assess with an unprecedented spatial resolution what is happening in vivo. This contributes to a better understanding of a myocardial perfusion at a microvascular level. In recent years, the development in this topic have led to a reduction of: a scan time, motion artifacts, use of a contrast agent and a radiation dose exposure. This has also led to a significant improvement in a spatial and a temporal resolution. However, all the above aspects still remain problematic and have their limitations.
For MPI, as mentioned above, we face the problem of artefacts that impede an accurate evaluation of a myocardial perfusion, namely artefacts that are caused by the motion of a heart. The artifacts may lead to the appearance of hypoperfused areas that may mistakenly resemble myocardial ischaemia. In that scenario, we are facing the problem of identifying what defects are relevant for perfusion and which are false positive findings. The above-mentioned motion artefacts can arise for several reasons and, for example, by a patient movement during the scan procedure or by a cardiac motion because of high or irregular heart rates. Examination of different cardiac phases is a way to reduce the impact of motion artefacts in MPI scans. Specifically, real (or true) perfusion defects can be observed in all phases of the cardiac cycle, whereas a motion artefact is visible only in a few cardiac phases. In addition, problems in a temporal resolution can be caused by the occurrence of motion artifacts, especially in patients with a high heart rate and a high heart rate variability. A sufficiently high temporal resolution as well as a sufficiently high spatial resolution are required during a dynamic CTP imaging of a heart rate at stress to image the entire myocardium which is “working” at a high rate (following Seitun S (2016) Stress computed tomography myocardial perfusion imaging: a new topic in cardiology, Revista Española de Cardiología (English Edition), 69 (2), p. 188-200, Seitun S (2018) CT myocardial perfusion imaging: a new frontier in cardiac imaging, BioMed research international, Volume 2018, Article ID 7295460, 21 pages). Cardiac patients, due to their medical conditions, may undergo multiple imaging procedures which increases their radiation exposure above what is expected for a lifetime. For example, the dynamic CTP is associated with high radiation doses, because this imaging technique requires multiple acquisitions to generate the TACs. During the static CTP, an average radiation patient exposure is in the range of 1.9-15.7 mSv with an average value of 5.93 mSv. For the dynamic CTP imaging, effective dose values are in the range of 3.8-12.8 mSv with an average radiation dose of 9.23 mSv. A CT-based dynamic myocardial perfusion imaging requires a relatively high radiation exposure. Various techniques are used to reduce radiation doses, for example, a low tube voltage is set during the dynamic CTP for patients with normal body mass indexes. In this way, the radiation dose can be reduced by 40% while maintaining image quality and a good MBF assessment. However, improvement of CTP imaging procedures is still needed. The improvement is to eliminate unnecessary radiation exposure and risks associated with this procedure (see Danad I (2016) Static and dynamic assessment of myocardial perfusion by computed tomography, European Heart Journal-Cardiovascular Imaging, 17(8), p. 836-844).
Methods based on the indicator-dilution flow analysis have been developed for a long time and are not only limited to application in typical medicine. There is a very wide range of applications of these method in hydrology and hydrogeology as well as in chemical and process engineering. In particular, the methods were used in cases of moving rocks or debris in a stream bed, when a stream cross-section was indeterminant or flow was physically inaccessible or unsafe (see Florkowski T, Davis T G (1969) The measurement of high discharges in turbulent rivers using tritium tracer. Journal of Hydrology, 8(3), p. 249-264, Kilpatrick F A, Cobb E D (1982) Measurement of discharge using tracers, U.S. Department of the Interior, U.S. Geological Survey, Evans G V (1983) Tracer techniques in hydrology, The International Journal of Applied Radiation and Isotopes, 34 (1), p. 451-475). One should note a change of terminology between engineering and medicine. In particular, two equivalent engineering terms “dye-dilution” and “tracer-dilution” are replaced by one medical term “indicator-dilution”. Accordingly, a cardiac output is measured in physiology and a flow discharge is measured in hydrotechnics. The latter was described by Charles Allen in 1927. Charles Allen named the method as the “salt-velocity” method (see Allen CM (1927) Hydraulic-turbine tests by the Allen method, Power Plant Engineering, 31(10), p. 549-551). Tetrach of Trachonitis (A.D. 37) used chaff to trace underground streams of the likely source of the Jordan (see page 247 in Assaad F A, LaMoreaux P E, Hughes T (2012) Field methods for geologists and hydrogeologists, Springer-Verlag Berlin and Heidelberg GmbH & Co. KG). In accordance with circulatory physiology, a mean transit time (MMT or an equivalently tissue transit time—TTT) is introduced. This results from the fact that a highly important parameter is residence time or residence time distribution in chemical engineering and reactor design (following Wingard L B, et al. (1972) Concepts of residence time distribution applied to the indicator-dilution method, J Appl Physiol, 33(2), p. 264-275). Nowadays, the methods based on indicator-dilution flow analysis (also classified as indicator-dilution class methods) are partially or fully substituted—mostly—by computational fluid dynamics based on numerical simulations (following Chen K (2018) CFD simulation of particle residence time distribution in industrial scale horizontal fluidized bed, Powder Technology, 345(1), p. 129-139, Walker P Sheikholeslami R (2003) Assessment of the effect of velocity and residence time in CaSO4 precipitating flow reaction, Chemical Engineering Science, 58(16), p. 3807-3816).
This review of prior art developments is underlying that there is a need for a reliable method that will allow to diagnose patients without the above-mentioned risks and problems. The present invention addresses this need.
The forthcoming description is to better understand the principle and benefits of the invention claimed in the appended claims. Therefore, it is not meant to be limiting in any sense. In particular, certain explanations to the meaning or to the function of, for example, equations used within the invention cannot be understood as limiting or exhaustive. They are provided to facilitate understanding and cannot replace the full understanding provided by a common general knowledge of a skilled person in this field.
The present invention redefines the myocardial perfusion imaging (MPI) paradigm. Its core and one of the biggest improvements over the paradigm is replacement of an in vivo imaging phase (usually realized with PET, SPECT, or MR scanners) with an in silico “imaging” phase while providing medically relevant and consistent results. The central element of the invention is a model of an arterial and myocardial blood flow solved using numerical fluid dynamic techniques. This way, all the insufficiencies of the known in vivo approach, including all the risks to the overall health of patients undergoing the procedure, are eliminated. In addition to that, a practically unlimited spatial-temporal resolution is achieved. Moreover, the result of the invention is completely free of motion artifacts and the patient is not additionally exposed to radiation. At the same time, the cost of performing this kind of in silico analysis is significantly lower.
Within the scope of this disclosure, the term “tracker” can be used interchangeably with the term “contrast”. The term “porous” should be understood as including or being equal to the term “saturated porous”.
In one aspect, the invention pertains to a method for assessment of a myocardial perfusion using a non-invasive measurement and poromechnics for a human patient. In particular, poromechnics involve computational fluid dynamics. The method can be realized as a computer-implemented method.
The method comprises generating an anatomical geometry model (
The anatomical geometry model (
In addition to above-mentioned benefits and the new way for approaching myocardial perfusion imaging (MPI), the use of the fluid domain and the porous domain allows for simplification of description of a myocardium connected with a plurality of coronary vessels. This allows for faster calculations and for arriving at perfusion metrics faster. At the same time, it should be emphasized that building an accurate complex model of a myocardium connected with a plurality of coronary vessels without the use of the fluid and porous domain may not always be possible.
In embodiments, the patient-specific anatomical structure data are obtained via at last one of a computed tomography angiography (CTA), a nuclear magnetic resonance (NMR) and/or a computed angiographic. This anatomical geometry model can be used in computational fluid dynamics to calculate perfusion metrics.
In an embodiment, the patient-specific anatomical structure data comprises images stored in a DICOM format. This allows for reconstruction of time attenuation curve (TAC) and arterial input function (AIF) curves for a number of time points.
The step of generating of the anatomical geometry model (
The method further comprises determining at least one patient-specific boundary condition and selecting at least one physical parameter for describing blood as a flow medium using patient-specific demographic and general medical data. The boundary condition and the physical parameter are to be used in computational fluid dynamics to, for example, calculate perfusion metrics. The term “patient-specific” means specific for patient and not e.g., data or parameters for a population of patients. The use of patient-specific data allows for obtaining results that are medically relevant for a particular patient and not mean results for a population of patients.
In one embodiment, the patient-specific demographic and general medical data comprise at least one of patient's age, sex, height, weight, systemic blood pressure, blood test results and/or information on the current patient's pharmacological therapy. The current patient's pharmacological therapy may be selected from, but not limited to, a beta-adrenergic blocking agent, an angiotensin-converting-enzyme inhibitor and/or an antiarrhythmic agent. In an embodiment, the physical parameter used to for describe blood as a flow medium is selected from the group comprising density and a dynamic viscosity coefficient.
The method further comprises simulating a blood flow in the anatomical geometry model (
In preferred embodiments of the invention, the fluid domain (Ωf) may satisfy the following equation
for which: ρ—fluid density, t—time coordinate, u and p flow velocity and pressure, η is dynamic viscosity coefficient.
In preferred embodiments of the invention, the porous domain (Ωp) may satisfy the following equation
for which: ρ—fluid density, u and p—flow velocity and pressure, η—dynamic viscosity coefficient, CF—inertial resistance tensor, K—permeability tensor.
In preferred embodiments of the invention, the flow through the fluid domain (Ωf) and the porous domain (Ωp) may satisfy the following equation:
for which χ=1 for the porous domain (Ωp), or 0 elsewhere,
For certain embodiments, in the porous domain, the Navier-Stokes equation are supplemented with additional terms associated with the resistance induced by a porous medium. The additional terms are based on the Forchheimer model describing both the viscous and inertial drag effects. This combination of the Navier-Stokes equation with the Forchheimer model allows to by-pass the extreme geometrical complexity of the porous domain, while it is able to recreate the nature of the phenomena taking place in said domain. In this inventive approach, information of key importance for diagnostic purposes, i.e., information about a flow distribution, is directly extracted. In other words, this unique feature allows to arrive at a relevant information without the time and effort needed to build and to solve a mathematical model of a porous domain. This is beneficial in terms of a need for computational power and in terms of pace of calculations.
The present invention allows for calculating of perfusion metrics which can include in embodiments:
The invention allows for calculating of MBF, MBV, MTT and TTP directly or in the form of indexes, ratios, regional values or patient-specific reference values.
Said perfusion metrics can be calculated using regional tracer kinetic recovery methods or using sequential Monte Carlo streamline trajectory Lagrangian tracking. In particular, calculation of the perfusion metrics can involve MaXimum Slope (MXS), Deconvolution of Impulse Response Function (DRF), multi-pathway, multi-species and non-linear Blood-Tissue eXchange (BTX). In other embodiment, calculation of the perfusion metrics can involve
wherein si are train vectors {Δs1, Δs2, . . . , Δsn}, where Δsi=f(t, ri).
Further characteristics and advantages of the invention will be more apparent from the detailed description of non-limiting embodiments and from the appended figures. It should be understood that embodiments presented in the description and in the claims, unless stated otherwise, can be combined together in any order and number to produce new embodiments that form part of the disclosure. The description comprises a great number of references to prior art documents, in particular to scientific journals. The full disclosure of said references is part of this disclosure.
The invention will now be described in more detail with reference to the appended figures in which:
While this description often refers to engineering as a source of medical approach used in the invention, there are significant differences between the engineering approach and the medical approach. Said differences are visible for a person skilled in the art and are not limited to different naming used in both fields. For example, certain engineering measurements can be made directly since, in a given example, we can measure a flow in a pipe or in a river and the shape of the pipe or the river does not change significantly during the measurement. On the other hand, a heart is undergoing cyclic significant changes but there is no indication in which part of a cycle a given measurement is made. This limited example is only to sketch the differences between engineering and medicine. This description of the invention is not meant to explain all aspects of engineering and medical convention used herein as this should be understandable for the skilled person.
The indicator-dilution methods used in cardiovascular diagnostics, as well as the well-known engineering dye-dilution, tracer-dilution or even the Allen's salt-velocity method, derive from one very elementary concept binding the three conventional scales: time, volume and mass. Labelling them—conventionally—as T, V, M we have T=V/Q (Q—volumetric flow rate), M/V=C (C—concentration of indicator) while in view of the change of M as a function of time we have the relation dM/dt=Q·ΔC (following Peters A M (1993) A unified approach to quantification by kinetic analysis in nuclear medicine. J Nucl Med. 34(4):706-13). The whole analysis concerning these three conventional scales and determined diagnostic parameters such as MBF, MBV, MTT, TPP by the software accompanying CT scanners is implemented according to equations (1) . . . (25). The result of an in vivo imagining comprises a series of images, usually saved in a DICOM format. The results and its particular format allow for reconstruction of TAC and AIF curves for a number of time points. Each time point is associated with dozens of images that must be saved by the device in order to reconstruct the tissue time attenuation curve. This can be technically demanding to be done within a split of a second or even within the time needed to record a single heart cycle. In fact, the time points on the timeline are quite remote from each other and only their mutual time interval is measured by tracking the patient's ECG (electrocardiogram). A representative example of this is
The idea of estimating the degree of circulatory insufficiency on the basis of an MPI examination is laudable and over time will certainly become a new gold standard for diagnostic tests. However, its implementation which is based on what was developed in engineering (e.g., in civil engineering or groundwater hydraulics) does not seem to lead to success. This can be explained. A heart is an organ that is subject to a constant motion and which is experiencing cyclic contractions. It is not possible, even for a few minutes, to stop a human heart for the needs of a diagnostic procedure. On the other hand, technical limitations of the diagnostic procedure do not allow to increase a time resolution significantly. Already at a very early stage of application of perfusion diagnostic methods, this problem was recognized. For example, a team from the University of Münster concluded: “Sampling intervals of >1 second on PCT imaging calculated using software relying on the maximum slope model significantly alter absolute CBF and TTP value . . . ” (see Kloska S P et al (2010) Increasing sampling interval in cerebral perfusion CT: limitation for the maximum slope model, Acad Radiol, 17(1), p. 61-66). Indeed, due to the interrelationship of time scales of a cardiac cycle and acquisition, and furthermore due to the magnitude of a volume of distribution, the derivative of dTAC/dt must be calculated from a section of the TAC curve which in turn is built using two, three, or four time point-layers. Even worse, the injection site for a contrast agent is usually far away from the target organ and the duration of the injection must be sufficiently long relative to the other time scales (which is required for the well-being of a patient).
In an embodiment of the present invention, it is permitted to use the maximum slope method (MXS) and, particularly, the relations (22)-(25) for recovery of regional tracer kinetics (TKR).
In another and preferred embodiment of the present invention, the recovery of regional tracer kinetics (TKR) will be based on an assessment of how the tracer is retained in the vasculature which is done instead of tracer rates. If, within a time interval Δt, a Δm(t) of a contrast agent from the total amount m introduced at the inlet has flows into the outlet region of interest (ROI), the empirical probability of such an event is expressed by
In this way, we can construct an empirical probability density function h(t) which is relating the probability value P(t∈Δt) with the width of the interval Δt
Knowing the distribution of the probability density function h(t) over the time over which fractions of the tracer substance reach the outlet, the empirical cumulative distribution function can be easily determined using the following equation
and the probability of the opposite event consisting of an event where time exceeds the assumed value of t
Mathematical description can be presented as: If a tracer in an amount m is introduced at the inlet in an infinitely short time interval, then at any time t an amount H(t)·m reaches the ROI volume, and there is still R(t)·m of the tracer remaining outside of the ROI. In the real-world conditions, the tracer is introduced in fixed portions, over a finite time interval Δti, often at a distal location, and it slowly flows with the blood to the presumed ROI. Then, the concentration change recorded by the TAC(t) curve induced by an earlier injection at the inlet at time Δti (which corresponds to the inlet value AIF(ti)), will be expressed by: TAC(t)=AIF(t)·T(t−ti)·MBF·ρt·Δti. If we treat the gradual changes in the concentration recorded at the arterial input function as a series of injections at times t1, t2, . . . , tn with durations Δt1, Δt2, . . . , Δtn, then at time t the concentration in the ROI volume will be expressed by
and consequently, with sufficiently large n→∞, the concentration in the target region of interest will change to
By introducing the impulse response function IRF(t)=ρt·MBF·R(t), it is easy to see that the resulting product integral defines a convolution operation which can be expressed by
The resulting equation is the central element of the secondary and preferred embodiment of the present invention. This equation defines a concise transformation principle of the state observed at the input to the state in the arbitrary region of interest and this may be expressed by:
Awareness of impulse response function IRF(t) greatly facilitates the assessment of perfusion parameters (perfusion metrics). Note that since max(R(t))=1, then according to the definition max(IRF(t))=ρt·MBF, hence the myocardial blood flow MBF=max(IRF(t))/ρt. By summing the equation (31) over a long period of time, by virtue of Fubini's theorem we arrive at the following equation
since MTT=∫0∞t·h(t)dt=∫0∞(1−H(t))dt=∫0∞R(t)dt, but MBF=MBV/MTT. The above-presented approach is part of DRF.
In a third preferred embodiment of the present invention, recovery of regional tracer kinetics (TKR) is carried out through a blood-tissue eXchange model (BTX) (following Wirestam R (2012) Using contrast agents to obtain maps of regional perfusion and capillary wall permeability, Imaging in Medicine, 4(4), p. 423-438, Bindschadler M et al. (2014) Comparison of blood flow models and acquisitions for quantitative myocardial perfusion estimation from dynamic CT, Phys Med Biol, 59(7), p. 1533-56, Tofts P S et al. (1999) Estimating kinetic parameters from dynamic contrast-enhanced T(1)-weighted MRI of a diffusable tracer: standardized quantities and symbols, J Magn Reson Imaging, 10(3), p. 223-32). The previously presented TKR models, in the context of pharmacokinetic and pharmacodynamic (PK/PD) modeling, are one-compartment models. Indeed, it would be easy to build a corresponding model of blood-tissue eXchange processes on top of a n-compartment PK/PD model. For example, a two-compartment eXchange model may be defined as follow
where c is tracer concentration (p—plasma, e—extracellular, extravascular space, ap—arterial plasma) v is fractional plasma volume in a compartment, PS is permeability-surface area product. When the effect of the tracer in the extravascular space is negligible, in an embodiment ct=vece and dct/dt=vedce/dt.
The present invention comprises in silico determination of perfusion parameters (perfusion metrics) MBF, MBV, MTT and TTP with the use of computational fluid dynamics (CFD) methodology. In a particular embodiment of the present invention, said determination is realized with the use of one of the regional tracer kinetic recovery (TKR) methods: MaXimum Slope (MXS), Deconvolution of Impulse Response Function (DRF), multi-pathway, multi-species, non-linear Blood-Tissue eXchange (BTX). The use of any of the listed methods is connected with an assumption of familiarity in time evolution of a marker concentration in the volume ROI in every point (or voxel), providing metrics TAC(t) being the result for arterial input. In the case of in vivo testing, time evolution of a marker concentration in the volume ROI in every point (or voxel) is provided straightforwardly by imaging results and presented in Hounsfield units (CFD does not present the same and CFD provides pressure and velocity distribution). After images co-registration and motion artifacts removing, a tissue attenuation curve in the volume of interest (VOI) voxel per voxel is reconstructed. This commonly adopted strategy is used by every producent of a perfusion diagnostic device because it is the easiest. The easiest strategy cannot cover all relevant cases and conditions. As it was mentioned earlier, factors that affect measurements are causing, especially in this case of this commonly adopted strategy, differences between the obtained results with different devices and/or with different protocols and/or with different cutoff values for MBF, MBV, MTT or TTP. The variation between the results provided by different researchers may by of dozens of precents (sometimes it is more than 100%). Additionally, the results obtained using this way do not have any direct physical interpretation which may be expressed in terms of e.g., established scales of flow, volume or time. Moreover, the results have really a loose connection with flow, volume or time as such. In consequence, more and more researchers are interested in the subject of in vivo perfusion imaging via advection-diffusion (see Liu P et al. (2021) Perfusion imaging: An advection diffusion Approach, IEEE Transactions on Medical Imaging, 40(12), p. 3424-3435). If with the flow rate any substance is transported, then its local concentration change ∂c(t, r)/∂t is described by advection-diffusion transport equation (Bird R B et al (2006) Transport phenomena, John Wiley & Sons) which may be expressed as
where q is source term, while flux (to be precise, a flux density as a rate of quantity per unit area) is expressed as
In a preferred embodiment of the present invention, marker concentration field distributions can be determined based on the above-mentioned advection-diffusion transport equation. To be more precise, by solving the equation numerically. In the case of an in vivo imaging results analysis, which is opposite to an in silico method being subject of the present invention, a velocity field does not have to be calculated but can be relatively easily reconstructed based on the imaging results. It is hard to estimate how this approach is going to be recognized by the producents of perfusion in vivo imaging devices. In any case, this approach can be very inspiring for development of in silico methods. In this case, a velocity field cannot be reconstructed but has to be calculated using a separate model. From the standpoint of engineering, the description of a marker flow in a tissue can be understood as a description of a marker flow in a porous medium or saturated porous medium. Usually, the description of this flow can implement a Henry Darcy equation (see Ingham D, Pop I (2005) Transport phenomena in porous media, Vafai K (2015) Handbook of porous media, CRC Press) which has a form of
where K is a permeability tensor. Considering that a velocity field satisfies a Euler balance equation (for liquids), a differential equation can be obtained in the form of
This equation gives the possibility of calculation a pressure distribution in a porous domain. Once we arrive at pressure (p), we can use the Darcy equation to obtain the velocity u (see Alves J R et al (2019) Simulation of the perfusion of contrast agent used in cardiac magnetic resonance: a step toward non-invasive cardiac perfusion quantification, Front Physiol, 10, p. 177, Alves J R et al. (2016) Simulation of cardiac perfusion by contrast in the myocardium using a formulation of flow in porous media, J. Comput. Appl. Math, 295(C), p. 13-24). Every differential equation case, including equation (36), requires boundary condition (here: pressure) to be solved. In every case—also in a preferred embodiment of the present invention—a pressure distribution on the edge of a porous heart muscle is affected by the inflow from coronary arteries. This observation can be generalized for perfusion models of different organs. In a simplified approach, one might use a Daniel Bernoulli equation to achieve this goal. (following Vladimír Lukes et al. (2014) Numerical simulation of liver perfusion: from CT scans to FE model, CoRR, abs/1412.6412). A more advanced approach to the modeling of an inflow to myocardium is presented in the work of a French-American team (see Papamanolis L et al. (2021) Myocardial perfusion simulation for coronary artery disease: a coupled patient-specific multiscale model, Ann Biomed Eng, 49(5), p. 1432-1447, Jaquet C (2018) Toward myocardium perfusion from x-ray computed tomography, Ph.D. Thesis, Université Paris Est) which developed the following equations
Said equations are describing the coronary flow as a steady blood flow along the centerline axis of each vessel. In fact, models implementing said equations are well known to any engineer studying available textbooks about hydraulics (see Wylie E B, Streeter V L (1967) Hydraulics Transients. McGraw-Hill). It is worth to note that this approach to the description of a coronary flow can be easily and effectively expanded by introducing inertial effects. A good example is the work of a Siemens Corporate Research and Technology team (see Sharma P et al. (2012) A patient-specific reduced-order model for coronary circulation, 2012 9th IEEE International Symposium on Biomedical Imaging, p. 832-835, Formaggia L et al. (2003) One-dimensional models for blood flow in arteries, Journal of Engineering Mathematics, 47, p. 251-276). Nevertheless, in its original form, so far, it has not been used in the context considered here. Equation 38 is not used in the invention.
In a preferred embodiment of the present invention, evolution of velocity (u) and pressure (p) fields is described by a system of Euler mass balance and Navier-Stokes momentum equations. A myocardium (and in general a muscle tissue) supported by a network of branched arteries can be considered to be two non-overlapping and simply connected domains
It is a common general knowledge that blood forms a polydisperse system comprising a liquid phase and a suspended dispersed phase. The liquid phase is plasma and the suspended dispersed phase is created from morphotic components: erythrocytes (which are the most numerous and responsible for a gas eXchange), leukocytes (whose numbers are low and which provide immunity functions) and thrombocytes (which are actively involved in coagulation processes, ensuring haemostasis). In average conditions, a blood density is in the range of ρf=1.05 . . . 1.06 g/cm3. A blood density is generally decreasing with age (reaching its maximum for newborns). A blood density can be more precisely determined as a function of a volume fraction of morphotic elements present in human blood. For example, the Hawksley formula is often used to estimate blood density: ρf=ρpl(1−HCT)+HCTρRBC, where ρpl, ρRBC is density of plasma and erythrocytes, and HCT is a hematocrit (see De Gruttola S et al. (2005) Computational simulation of a non-newtonian model of the blood separation process, Artif Organs, 29(12), p. 949-59). Estimation of a blood viscosity coefficient is much more complex task in comparison to a blood density (following Milnor W R (1982) Hemodynamics, Williams & Wilkins, Vlachopoulos Ch et al. (2011) McDonald's blood flow in arteries, CRC Press, Kannojiya V et al. (2021) Simulation of Blood as Fluid: A Review from Rheological Aspects, IEEE Rev Biomed Eng, 14, p. 327-341). First of all, the coefficient depends on the content of the solid fraction, which is considered to be typical for dispersion systems. For blood, the coefficient is expressed as the value of the hematocrit number (HCT). Moreover, red blood cells have a unique ability to create three-dimensional structures in the form of rouleaux (see Samsel R W, Perelson A S (1982) Kinetics of rouleau formation. I. A mass action approach with geometric features, Biophys J, 37(2), p. 493-514, Samsel R W, Perelson A S (1984) Kinetics of rouleau formation. II. Reversible reactions, Biophys J, 45(4), p. 805-2). In a preferred embodiment of the present invention, blood rheological properties to be used in computational fluid dynamics are described using a model which is of the Williamson class (Williamson fluid model), for which (η−η∞)/(η0−η∞)=F(HCT, {dot over (γ)}). The structural function of the F model in a preferred embodiment is selected according to the Carreau, Cross, Yasuda model (see Carreau P J (1972) Rheological equations from molecular network theories, Journal of Rheology, 16(1), p. 99-127, Cross M M (1965) Rheology of non-Newtonian fluids: A new flow equation for pseudoplastic system, Journal of Colloid Science, 20(5), p. 417-437, Yasuda K (1979) Investigation of the analogies between viscometric and linear viscoelastic properties of polystyrene fluids, Ph.D. Thesis, Massachusetts Institute of Technology). In another embodiment, it is also possible to use the Walburn-Schneck model, which explicitly takes into the account the influence of blood plasma proteins (see Walburn F J, Schneck D (1976) A constitutive equation for whole human blood, Biorheology, 1, p. 201-210).
While it is customary for engineers to express the movement in porous domain Ωp by the Darcy equation (see eqn. (36)), this equation has been specifically rejected and it is not implemented in the present invention. This equation does not include inertial effects in its structure (in principle, it can only be used within the Reynolds Rep<1). Instead, a non-linear equation is used in a preferred embodiment of the present invention (following Forchheimer P (1901) Wasserbewegung durch Boden, Z. Ver. Deutsch. Ing, 45(1901), p. 1782-1788, Whitaker S (1996) The Forchheimer equation: A theoretical development, Transport in Porous Media, 25, p. 2761, Vafai K (2015) Handbook of porous media, CRC Press) which has the form of
where CF is an inertial resistance tensor. In order to obtain a consistent phenomenological description, in a preferred embodiment of the present invention penalization approach has been introduced to model the flow through parts of the fluid Ωf and porous domain Ωp (following Bruneau C H, Mortazavi I (2008) Numerical modelling and passive flow control using porous media. Computers & Fluids. 37(5):488-498) with the following equation
while χ=1 in a porous domain, or 0 elsewhere.
Since the solution of the above equation provides the full description about the evolution of the pressure fields and about the velocity (in the nodes of mesh used for calculations) and thus about the local values of flow rates in each mesh cell, in another embodiment of the present invention, a myocardial blood flow (MBF) can be calculated directly. In another embodiment of the present invention, the reconstruction of the parameters describing perfusion can be realized with the use of sequential Monte Carlo streamline trajectory Lagrangian tracking in the form of
by train vectors {Δs1, Δs2, . . . , Δsn}, where Δsi=f(t, ri) (see Crane M, Blunt M (1999) Streamline-based simulation of solute transport, Water Resour. Res, 35, p. 3061-3078, Salamon P et al. (2006) A review and numerical assessment of the random walk particle tracking method, J. Contam. Hydrol, 87, p. 277-305, Friman O et al. (2010) Probabilistic 4D blood flow mapping, International Conference on Medical Image Computing and Computer-Assisted Intervention, 13(Pt 3), p. 416-23, Bear J (2018) Modeling phenomena of flow and transport in porous media, Springer International). This form of post-processing of the calculation results enables to determine distributions of the probability density function h(t), the distribution function H(t) and the residual function R(t), in accordance with the equations (27)-(29). In this way, it is possible to use any of the above-mentioned methods for calculating perfusion parameters based on the results of CFD simulations.
An embodiment of the method of the present invention comprises six steps: 1) collection of patient-specific demographic and general medical data and patient-specific anatomical structure data, 2) quantification of boundary condition, 3) creation of a discrete model of geometry, 4) numerical simulation in fluid and porous domains, 5) recovery of regional perfusion metrics, and—finally—6) visualization of the regional perfusion metrics.
In the first step, patient-specific demographic and general medical data are collected. In particular, patient-specific demographic and general medical data include age, sex, height, weight, systemic blood pressure and/or blood test results (if available, hematocrit number in particular) of a particular patient. The patient's demographic and general medical data may further include information on a patient's pharmacological therapy. The patient's pharmacological therapy may include a beta-adrenergic blocking agent, an angiotensin-converting-enzyme inhibitor and/or an antiarrhythmic agent. Other embodiments comprise other pharmacological therapies. It should be noted that pharmacological therapies include also different dosing regimes as well as any further medical effects of drug therapies. It is contemplated that in certain embodiments only selected patient-specific demographic and general medical data are collected, e.g., gender, age or pharmacological therapies. The patient-specific demographic and general medical data can be collected via an interview with a patient or can be obtained from a database.
In addition to patient's demographic and general medical data, which are subsequently given a physical interpretation, patient-specific anatomical structure data are also collected. In one embodiment, the patient-specific anatomical structure data comprise the results of a computed tomography angiography (CTA) of the coronary arteries and myocardium. In other embodiments, computed tomography and/or magnetic resonance techniques can be used. The patient-specific anatomical structure data comprise information about anatomy of a heart myocardium and vessels connected to the myocardium. The patient-specific anatomical structure data can be collected before, after or in parallel with the collection of the patient-specific demographic and general medical data. The patient-specific anatomical structure data can be retrieved from a database.
In the second step, the patient's demographic and general medical data is used to determine one or more patient-specific boundary conditions, as well as to select one or more physical parameters that will be describing blood as a flow medium. This is needed for a numerical simulation. The selected physical parameters may include a density and/or a dynamic viscosity coefficient. In the diagram of
As an example,
The resulting surface model is discretized through an unstructured volumetric mesh which is shown in
In third step, a computational grid with the conditions set on its edge (boundary conditions) and with the specification of the fluid together with porous domain physical properties are transferred to a solver engine. As presented earlier, two sets of equations are solved, albeit in a penalized approach. The equations of fluid motion in the form of a system of Euler and Navier-Stokes equations are solved for the entire computational grid. A skilled person will understand that the entire computational grid excludes boundary conditions which are put to arrive at a specific solution. In a porous domain, the Navier-Stokes equation are supplemented with additional terms associated with the resistance induced by a porous medium. The additional terms are based on the Forchheimer model describing both the viscous and inertial drag effects. This combination of the Navier-Stokes equation with the Forchheimer model allows to by-pass the extreme geometrical complexity of the porous domain, while it is able to recreate the nature of the phenomena taking place in said domain. In this inventive approach, information of key importance for diagnostic purposes, i.e., information about a flow distribution, is directly extracted. In other words, this unique feature allows to arrive at a relevant information without the time and effort needed to build and to solve a mathematical model of a porous domain. This is beneficial in terms of a need for computational power and in terms of pace of calculations. At the same time, it should be emphasized that building an accurate complex model of a heart myocardium connected with a plurality of coronary vessels without the use of a fluid and porous domain may not always be possible.
In the fourth step, a diagnostic interpretation by algorithms is given to the results obtained with the use of numerical fluid dynamics methods. Note that, regardless of the type and degree of complexity, the numerical fluid dynamics models only provide information about the velocity and pressure field distributions in the flow. These results are not directly used in medical diagnosis. When assessing perfusion, the overall expectation is that the diagnostic method will provide some contractual metrics such as MBF, MBV, MTT, TTP, or somehow related thereto (such as relative index values). The fourth step can be carried out in three variants, differing in computational inputs and achieved results. In the most direct approach (
In fifth step, reconstruction of the tracer kinetics is performed for the determination of standard perfusion metrics and it is related to the results obtained in step 4.B or 4.C (results from 4.A are not subjected to this kind of processing). It is possible to use any of the tracer kinetic recovery (TKR) strategies discussed earlier in details, i.e., MaXimum Slope (MXS), Deconvolution of Impulse Response Function (DRF), Blood-Tissue eXchange (BTX). The final results are exported for later visualization.
In sixth and last step, only the visualization of previously obtained perfusion metrics is performed. As shown, the aforementioned
Any calculation step or substep of the method according to the invention can be implemented using a computer or a computer program. In some specific embodiments, some or all calculations are done using a computer program stored on a computer or on any type of a memory device or both. In another embodiments, some or all calculations for the purpose of the method can be done remotely e.g., using cloud-based infrastructure which can include the use of Internet or a local network.
While all measurement methods and algorithms described here are intended to define the parameters of the invention and to provide tangible results to be compared with clinical trials, there are by no means limiting and are exemplary. The words like “including” or “comprising” are not limiting and, for example, when an element A includes another element B, the element A may include other element or elements in addition to the element B. The use of a singular or plural form is not limiting for the scope of the disclosure and, for example, a part of the description which is indicating that an element A contains an element B, this part also discloses that multiple elements B are contained in one element A and multiple elements A are contained in one element B as well as an embodiment where multiple elements A contain multiple elements B. Many other embodiments will be apparent and clear for those skilled in the art upon reviewing the content of the present disclosure. The scope of the invention, therefore, should be determined with reference to the appended claims along with the full scope of equivalents to which said claims are entitled to.
| Filing Document | Filing Date | Country | Kind |
|---|---|---|---|
| PCT/IB2022/051883 | 3/3/2022 | WO |
