本申请要求2010年6月24日提交的美国临时申请系列号61/358,400的权益,且也含有可能与在下述申请中公开的和/或要求保护的主题有关的主题:2007年3月22日提交的美国专利申请号11/575,846、2007年3月27日提交的美国专利申请号11/691,748、2007年3月27日提交的美国专利申请号11/576,060(现在授权为美国专利7,925,330)、2009年6月12日提交的美国专利申请号12/519,040、2009年6月15日提交的美国专利申请号12/519.213、和2010年1月15日提交的美国专利申请号12/669,292,它们的公开内容通过引用并入本文,并形成本文的一部分。
具体实施方式
如在本文中和在所附权利要求书中使用的,单数形式“一个”、“一种”和“所述”包括复数指示物,除非上下文另外清楚地指明。因而,例如,对“一个参数”的提及包括多个这样的参数和本领域技术人员已知的它们的等效参数,诸如此类,对“所述参数”的提及是对一个或多个这样的参数和本领域技术人员已知的它们的等效参数的提及,诸如此类。
在几个代表性的实施例中,使用生理学基础上药代动力学或PBPK建模范例。在许多实施例中,所述PBPK模型可以例如转化成离散时域。此外,所述模型的许多实施例会提供运输延迟在药物动力学中的明确导入(也就是说,传输时间,例如,造影剂穿过肺脉管系统的传输时间)。可以描述药物在遍布于体内的血管结构中的分布。尽管本文描述的模型可以用于模型化药物的分布,通常,在几个代表性的实施例中,所述模型的主题是,在对比分布的动脉相过程中的CT血管造影术。因为本文描述的模型组合了许多建模拓扑学,且包括Bae对比模型的方式的参数化,这样的模型在本文中统称为“混合”模型。但是,如上所述,不同于例如Bae模型,所述混合模型可以转化成离散时域。此外,除了将可构建的运输延迟术语整合进所述模型结构和输入(外周静脉的)隔室中的饱和非线性内以外,所述混合模型的许多实施例会提供不同造影剂浓度的注射相模拟,所述注射相包括在其中注射几乎不含有或不含有对比增强剂的流体的一个或多个注射相(例如,稀释剂或盐水注射相,其有时称作在造影剂注射以后的冲洗或驱逐相)。以前的(例如造影剂传输的)生理学基础上的药代动力学模型,不会模型化例如这样的盐水或其它“冲洗”相(其通常用于例如CTA成像操作过程中)对体内对比增强的影响。本文中的混合模型能够实现例如造影剂运输的模型化,甚至在造影剂注射已经结束以后,因为造影剂或盐水或不包含对比增强剂的其它快速推注的体积流率会增加周边隔室中的血液的体积流率。独立地处理血液的体积流率和对比增强剂的质量流率。在许多实施例中,所述混合模型包括:与例如整体Bae模型相比,减少的状态数目,而混合模型和整体Bae模型的输出之间的对比是有利的。例如,通过使用混合模型,测得更低的预测误差。当使用临床数据来对比2种模型的性能时,在某些情况下,所述混合模型胜过Bae模型。
在许多实施例中,使用数据驱动的方法来预测或模型化药物传输,其中使用一条或多条对比时间增强曲线(TEC),所述曲线得自例如造影剂的测试注射。在许多实施例中,使用基于模型的或参数识别技术来鉴别PBPK模型(诸如上面讨论的混合模型)的一个或多个参数。在许多这样的实施例中,采用模型简化策略,其中在造影剂的首次通过(例如,外周静脉室、右心隔室、肺隔室和左心/主动脉隔室)过程中考虑所述混合模型的最相关的隔室。在本文的代表性的研究中,考虑到可从测试快速推注TEC得到的数据的有限持续时间,没有尝试拟合所述混合模型的所有参数。但是,在提供足够数据的情况下,可以鉴别所述混合模型的所有参数。
还可以使用非参数的或不依赖于模型的鉴别技术来产生药物传输的非参数估测值。在几个实施例中,通过使用截短的奇异值分解(tSVD)解卷积技术求解反转的问题,产生药物/患者系统的非参数估测值。
本文所述的模型不仅可以用于预测药物(诸如造影剂或物质)在人心血管系统中的药代动力学,而且可以提供用于产生或计算注射操作的一个或多个参数(例如,注射方案和/或扫描的一个或多个参数)的系统和/或方法,所述注射操作会实现例如单个患者和操作的预期中的选择的增强靶标,同时使用减少的或最小体积的造影剂。
例如,在许多实施例中,如上所述的数据驱动的对比增强预测方法被用于方案/参数产生方法、系统或算法中,以成功地产生在一系列操作和患者变量之间的对比增强测线。
如本文关于注射操作使用的,术语“注射方案”或“方案”表示一组注射变量或参数,诸如流率、注射体积、注射持续时间、造影剂浓度等,它们限定了例如在注射操作过程中要递送给患者的流体的时机、量和/或性质。这样的参数可以在注射操作进程内变化。本文使用的术语“相”通常是指一组参数,它们限定了例如在小于注射操作的总持续时间的时间段(或相持续时间)内要递送给患者的流体的时机、量和/或性质。因而,相的参数会提供与所述相的持续时间相对应的时间实例内的注射描述。特定注射操作的注射方案可以例如描述为单相(单个相)、二相(2个相)或多相(2个或更多个相,但是通常超过2个相)。多相注射也包括这样的注射:其中所述参数可以在注射操作的至少一部分内连续变化。
可以确定的扫描仪参数包括、但不限于:向患者传送的辐射的量,能量输入(例如,电压或电流)、定时(例如,扫描开始时间、停止时间、延迟时间和/或持续时间)。
在几个实施例中,包含本文所述模型的系统可以包括这样的注射系统,其包括、例如,在图1B中描绘的双注射器注射器系统100。双注射器注射器系统公开在,例如,美国专利号6,643,537、公开的美国专利申请公开号2004-0064041和PCT国际专利申请号PCT/US2007/026194)中。注射器系统100可以例如包括双流体递送源(在本文中有时称作源“A”和源“B”;诸如注射器),其可运行以将第一流体和/或第二流体(例如,对比增强流体、盐水等)独立地(例如,同时地,同时地以彼此不同的体积流量比例,或顺序地,或彼此相继地(也就是说,先A再B,或先B再A))引入患者中。在图1B的实施例中,源A与增压机构(诸如驱动部件110A)可运行地相连,源B与增压机构(诸如驱动部件110B)可运行地相连。所述注射系统包括与注射器系统100可运行地相连的控制系统200,所述控制系统200可运行以控制驱动部件110A和110B的运转,从而分别控制来自源A的流体A(例如,造影剂)的注射和来自源B的流体B(例如,盐水)的注射。控制系统200可以例如包括用户界面或与用户界面通信,所述用户界面包括显示器210。在图1B描绘的实施例中,描绘了显示屏的一个实施例的一部分,所述显示屏显示了注射流率、注射体积和注射持续时间的参数区域,例如,关于流体A和/或流体B的3相注射。使用本文的参数产生系统和方法,可以产生一个或多个这样的相的参数。图1C阐明了显示屏的另一个实施例。
可以给用户提供选项,以调节和/或废除产生的方案或参数(例如,经由手动输入系统205,包括计算机领域已知的小键盘、键盘、鼠标等)。控制系统200可以包括与存储器或存储系统230可运行地相连的处理器220(例如,本领域已知的数字微处理器)。
本领域技术人员显而易见,许多注射器或流体递送系统,包括在例如美国专利号7,326,186、7,094,216、6,866,654、6,972,001、6,699,219、6,471,674、6,306,117、6,149,627、6,063,052、5,920,054、5,843,037、5,827,219、5,739,508和5,569,181中公开的多患者流体递送系统,也适合与本文的模型一起使用。
成像系统300可以例如是如上所述的CT系统、磁共振显像仪(MRI)系统、超声成像系统或正电子发射断层扫描术(PET)系统或单光子发射计算机体层摄影术(SPECT)系统。注射器系统可以与成像系统300通信连接,或与成像系统300部分地或完全地集成。成像系统300和注射器系统100可以例如经由本领域已知的输入/输出端口(由图2B中的箭头的终点表示)而通信连接。在图1B中,将成像系统300和注射器系统100例如解释为通过共同的通信集线器400而通信连接。或者,可以建立直接通信链路。使用计算机领域已知的一种或多种手动输入系统(例如,小键盘、键盘、鼠标等),可以手工地输入来自成像系统300和注射系统100之一的其它数据。还可以按照在例如公开的PCT国际专利申请号WO 2008/011401(其公开内容通过引用并入本文)中所述,部分地或完全地集成成像系统300和注射器系统或注射器100。或者,注射系统和成像系统300的所有解释的组件之一或一些,也可以与另一个单独的组件(其安置成与其它系统组件通信连接)集成或整合在后者内。
实现本文所述的系统和方法或其任意部分的软件和/或硬件可以例如被包括或整合在所述系统的一个或多个组件内(例如,在注射器系统100内和/或在成像系统300内),或在一个或多个由系统500表示的单独的或独立的系统内,所述系统500可以例如包括至少一个处理器(例如,数字微处理器)、存储系统520、显示器510和手动输入系统505。在图1B阐明的实施例中,系统500显示成与通信集线器400通信连接。如上所述,还可以建立直接通信链路。使用计算机领域已知的一种或多种手动输入系统(例如,小键盘、键盘、鼠标等),可以将来自一个或多个系统的其它数据手工地输入一个或多个其它系统中。实现本发明的系统和方法的软件(包括、例如,其一种或多种可执行的计算机算法)可以例如储存在存储器530中,并由处理器520执行。本领域技术人员显而易见,所述模型化和/或参数产生方法和/或系统的所有或部分功能可以替代性地存留在成像系统300(其可以例如包括至少一个处理器320、存储系统330、显示器310和手动输入系统305)中和/或在注射器系统100中。
在几个代表性的实施例中,本文所述的方法和系统与心胸脉管系统的CT结合使用。但是,本文关于患者特异性的心胸脉管系统的对比增强CT所述的方法和/或系统,也可应用于身体的其它解剖学区域。具体地,作为现代CT扫描仪的获取速度的结果,腿动脉的外周动脉血管造影术是挑战性的。在许多情况下,扫描仪必须减慢,以解释穿过脉管系统运输造影剂快速推注的生理学过程。通常,需要外周动脉CTA研究的患者也是糖尿病患者,或具有其它肾功能不全,且由于造影剂的大体积而潜在地易发肾损伤。因此,如本文所述计算出计算出所述患者中的减少的或最小的造影剂体积,可以是例如合乎需要的。
本文所述的方法和/或系统还可以用于神经学CT成像中。造影剂的简短快速推注(精确地且单个定时地),对于脑动脉的CTA而言是重要的。在这些研究中,重要的是,在造影剂填充头静脉之前,使扫描同步化。
本文所述的方法和/或系统还可以用于预测性对比增强情况中,其中不施用测试快速推注。有些辐射学家对造影剂跟踪软件的偏好胜过测试快速推注方法,以使在造影剂到达时扫描的获取同步化。在这些情况下,所述混合模型(假定患者人口统计数据是可得到的)可以用于先确定造影剂注射方案,以迭代方式达到希望的靶标。当构建对比方案时,患者的真实血流动力学状态通常是不可得到的。并且,因为造影剂在心肺回路中的传输的运输延迟,一旦注射造影剂,对特定患者的造影剂不再有“控制”。尽管如此,在该开环方案中使用所述模型仍然优于对患者的体型或生理学不做任何考虑的给药。此外,可能在模型预言性对照框架中使用所述混合模型和在扫描仪的造影剂跟踪获取过程中收集的数据。
在MRI检查过程中,常规地递送基于钆(Gd)的造影剂,以提供增强的动脉和静脉结构显影(MR血管造影术)。本文所述的方法和/或系统可以例如用于MR(和其它成像用途),同时例如认识到,(例如,Gd造影剂的)信号强度和血液-血浆浓度之间的关联是非线性的,正如CT造影剂一样。
混合模型
生理学基础上的药代动力学建模(PBPK)是一种建模方案,其在确定模型结构时会考虑有关的生理学和功能。在PBPK模型中,将身体分成许多互联的、与解剖学区域相对应的隔室。基于生物体的生理学和解剖学考虑,将每个隔室参数化(体积、血流量、灌注),其各自通过血管隔室相连,所述血管隔室会促进物质进出隔室的传递透过。该方案的一个优点是,在不改变模型结构的情况下在物种之间缩放所述模型的能力。PBPK模型的数学基础是各个隔室之间的质量守恒。
本文的造影剂混合模型的一个实施例的隔室显示在图1A中。所述模型中的每个子系统代表身体的一个解剖学区域。在PBPK模型中的子系统被分成3个隔室——细胞内空间、细胞外空间和血管内空间。造影剂不会进入细胞内空间,所以它在造影剂传输的PBPK模型的许多实施例中被忽略。
PBPK混合模型的几个实施例的状态变量x是隔室中的造影剂的质量(xi)。第i个隔室中的血液的体积流率/造影剂用Qi表示,体积用变量Vi表示。造影剂从子系统中的清除Cl,通过不可逆过程而发生。通过肾的肾小球过滤,发生造影剂的提取(Extraction)。
在其中意图使用模型来研究和预测造影剂在CT血管造影术过程中的分布的实施例中,定量造影剂在体内的首次通过过程中的血浆浓度是首要的问题,关于描述造影剂在扩散进全身器官和薄壁组织(除了肺子系统以外)中时的吸收和分布方面的兴趣更少。在这样的实施例中,所述模型仅需要考虑外周静脉的血管隔室、右心和左心以及全身循环。
[110]进或出模型隔室的质量“通量”(J)是:
在方程式1中,Q是体积流率,V是隔室的表观容积。流量系数fBODY,和fPER会缩放穿过各个身体段的心输出量。流量系数的总和必须等于心输出量QCO。fPERQCO+fBODYQCO=QCO。对于该模型中的血管隔室,所述表观容积是血管内的血液体积。将质量平衡应用于所述模型的子系统,会得到图1A中的子系统(外周、肺、右心、左心和身体)的下述总表达式(i=1:5)
在方程式2中,uexog(t)是造影剂的外源施用。它仅在周边隔室中是非零的,且定义如下:
uexog(t)≡Cinj(t)·Qinj(t) (3)
因为施用流率Qinj(t)或施用的造影剂的浓度Cinj(t)可以作为时间的函数而变化。体积流率和造影剂浓度的乘积称作碘施用速率。
将造影剂注射进人的左臂或右臂的外周静脉中。臂的外周静脉中的血液引流的内在流率范围是2-4ml/s。造影剂在小外周静脉中的高流率可能产生血流量,并因而,将注射速率加入外周静脉的血液的内在流率中。CTA检查的典型流率大于穿过外周静脉的血液引流的内在流率。将造影剂注射进外周静脉中,会导致穿过静脉的内在流率的累加贡献。在数学上,外周静脉子系统是外周子系统中的质量平衡的线性的、随时间变化的公式,那么,由于造影剂的注射:
with 其中
QPER(t)=QPER_END+Qinj(t) (4)
在方程式4中,QPER_END是穿过外周静脉的血液的内在流率,TBODY是在整个身体中再循环回外周静脉的造影剂的运输延迟。为了避免公式化整个模型的随时间变化的模型,可以重构方程式4,其中考虑到,在注射以后,QPER(t)=QPER_END,并且,在注射过程中,流率是QPER_END+QInj的总和。通过循环系统引入的运输延迟由TPER和TBODY表示,并代表来自外周和身体子系统的造影剂的输出延迟。时间延迟可以是常量,或是各个系统中的血液的体积和流率的函数。下面进一步讨论了关于时间延迟的细节。在随后的开发中,Tinj代表造影剂注射持续时间。在这些假设下,外周子系统动力学可以在标准的LTI状态空间公式中表示为:
在方程式5中,yPER(t)是离开外周子系统的造影剂的质量通量(单位gI/s)。外周子系统中的浓度是xPER(t)/外周静脉中的血液体积。
方程式4的通解是:
饱和行为会影响进入周边隔室中的注射流率。实验证据提示,大于8-10ml/s的注射会饱和,因为外周静脉的适应性质,或因为造影剂穿过右心房回流进下腔静脉中。因此,周边隔室动力学的一般描述是:
其中,非线性的外源输入函数是:
uexog(Qinj,t)=Qinj(t)·Cinj(t) Qinj(t)≤8[ml/s]
uexog(Qinj,t)=8·Cinj(t) Qinj(t)>8[ml/s] (8)
应用于图1中的右心子系统的质量守恒是:
肺子系统与纯血管子系统的差别在于,将毛细血管床和肺之间的渗透作用建模。为了完整性,且针对其中关心造影剂在组织中的再循环和积累作用的实施例,考虑组织隔室。
在方程式10中,QCO是心输出量,kLUNG_BT是造影剂从血液隔室向组织隔室中的速率转换系数,kLUNG_TB是从组织隔室返回血液隔室的速率转换系数,和CLLUNG_T代表离开组织隔室的不可逆清除期。
将穿过肺脉管系统的造影剂的延迟模型化为左心子系统的动力学中的输入延迟术语:
且yLUNG(t)是离开肺子系统的造影剂的质量通量。TLUNG是穿过肺脉管系统移动的造影剂的延迟时间。它可以是标量常数,或是心输出量和肺的血液体积的函数。
就CT血管造影术应用而言,造影剂快速推注的首次通过分布和传输具有重要意义。因为器官、肌肉和脂肪隔室会影响造影剂在再循环几次以后的分布,它们在几个实施例中未予考虑。其他人已经将全身循环组合为一个大的血管隔室,并以在许多实施例中的方式将全身循环建模。
可以阐述所述模型的统一公式(其组合了方程式4-12),以易于在数字模拟中实现。为了标记方便,并为了确保运输延迟仅作为输入和输出延迟而出现,定义了增强的状态向量,其包含所述模型的每个部分的状态变量。将右心、肺和左心子系统组合进心肺(CP)状态向量中:
xAug=[xPER xCP xBODY] (13)
其中所述心肺(CP)状态向量由右心、肺和左心状态变量组成:
xCP=[xRH xLUNG xLH] (14)
总系统表示为:
yPER(t)=CPERxPER(t)
yCP(t)=CCPxCP(t)
yCP(t)=CBODYxBODY (15)
外周子系统的输入向量由穿过外周静脉的外源性造影剂注射和来自身体子系统的再循环的造影剂组成:
流量系数fPER会缩放身体子系统中的造影剂的血液-血浆浓度。在方程式16中的标量函数φ(uexog)定义了在方程式8中描述的注射部位和右心之间的静脉系统的非线性饱和行为。
状态矩阵是:
其中穿过外周循环的随时间变化的流量QPER(t)由方程式4定义。对照矩阵(B’s)是:
BPER=[1 1]T (21)
BCP=[1 0 0 0]T (22)
BBODY=[1] (23)
所述模型的全积分仪视图的一个实施例显示在图2中。
在为数值模拟选择模型参数值时,可以在保真度和合理的可实现性之间做出平衡。难以推论出与个体的真实生理学参数匹配的、每个亚隔室的参数值。在许多实施例中,使用标准的生理学“查找”表和通过对群体数据的回归分析确定的关联来确定一个或多个参数。例如,可以使用任意数目的关联来确定一个或多个参数。如下面讨论的,还可以使用来自患者的测试快速推注增强数据的各种PBPK模型辨识方法来确定参数。在许多模拟研究中,使用仅仅基于人口统计数据(身高、体重、性别、年龄)的参数值和操作特异性的值诸如造影剂的流率、体积和浓度。
如在Guyton,A.,Circulatory physiology:cardiac output and itsregulation.1963,Philadelphia PA:Saunders中所述,针对年龄校正过的回归公式可以用于估测心输出量和血液体积。在许多实施例中,使用这里的心输出量估测,并在模拟所述模型时提供心输出量估测。所述估测器是:
其中参数h、w、和a是身高[英寸]、体重[磅]和年龄[岁]。在公开的美国专利申请号2010/0030073中,也讨论了心输出量的估测或测量。同样地,从公开的回归公式导出所述模型中的中枢血液体积估计值,且是身高、体重和性别的函数。为男性估测的血液体积是:
为女性估测的血液体积是:
在所述模型中的局部血液体积参数等于通过方程式28或方程式29估测的总血液体积。穿过每个亚隔室(除了外周以外)的血液的体积流率是通过方程式27计算出的心输出量。例如,根据在Bae,K.T.,J.P.Heiken,和J.A.Brink,Aortic and hepatic contrast medium enhancement at CT.Part I.Prediction with a computer model.Radiology,1998.207(3):第647-55页中所述的关联,可以设定在不同子系统中的血液体积。但是,不是象Bae等人那样将左心血液体积设定为总血液体积的3.6%,在本文的PBPK模型的几个实施例中加入额外的100ml血液,以在计算中包括升主动脉。在许多实施例中,用于构造给定受试者的血液体积的操作是,使用身高、体重、年龄和性别,计算心输出量和血液体积的估测值。然后使用表1中的关联,计算局部血液体积。
表1.在新颖的模型模拟中使用的局部血液体积参数
参数 |
值[单位] |
VPER |
.08*BV(身高,体重,性别)[ml] |
VRH |
.036*BV(身高,体重,性别)[ml] |
VLB |
.088*BV(身高,体重,性别)[ml] |
VLH |
.056*BV(身高,体重,性别)[ml] |
VSYS |
.848*BV(身高,体重,性别)[ml] |
其它模型参数包括肺、血液和组织隔室与中枢血液隔室的造影剂清除之间的速率转换系数。从血液供给和组织隔室提取造影剂发生在快速推注后数分钟,清除期(CLUNG_T)可以设定为零。同样地,在造影剂首次通过时,造影剂几乎没有穿过肺毛细血管床,且随后在组织中的积累是微小的。因此肺转移速率常数可以设定为零。
考虑的最后的参数是穿过外周静脉循环的造影剂快速推注传输延迟(TPER)、穿过肺系统的传输时间(TLUNG)和遍布整个循环系统的造影剂的再循环延迟(TBODY)。
在许多实施例中,系统再循环延迟保持恒定在例如30秒。但是,就其它2种传输延迟而言,考虑一个或多个全患者(per-patient)或患者特异性的变量或值。TPER是外周静脉子系统中的血液体积以及注射速率和内源静脉流量的总和的函数。同样地,穿过肺隔室的传输延迟是每位受试者的肺血液体积和估测心输出量的函数。还可以根据一个或多个患者特异性的变量来做出系统再循环延迟。在几个实施例中,为TPER和TLUNG延迟做出平推流的假设。在所述模型中使用的通过延迟的描述给出在表2中。
表2通过延迟参数
下面呈现了所述模型作为各种参数和输入的函数的性能,并与通过Bae模型做出的预测相对比(二者公开在且源自Simulink中的模型的实现)。描述了新模型的预测人数据集中的对比增强的能力,并也与Bae模型预测相对比。
在许多模拟研究中,改变注射参数和模型参数以证实所述模型的模仿已知的造影剂动力学行为的能力。进行模拟,所述模拟证实了所述模型的再现改变的心输出量对对比增强的影响的能力。通过模型模拟,证实了随着注射速率增加而发生的对比增强的饱和行为,正如造影剂注射以后的盐水冲洗的影响。
图3呈现了具有标称属性的虚拟患者的肺动脉和升主动脉中的模拟对比增强。实现所述模型,并在MATLAB(R2008b)中模拟。该实施例的患者参数是:体重=170磅,身高=68英寸,性别=男性,年龄=35岁。使用方程式27,估测的心输出量是6.60L/min。造影剂注射方案包括使用370mgI/ml浓度造影剂的对比,在5ml/s的体积流率(Qinj(t))注射20秒(体积=100ml)。在5ml/s,注射30ml体积的盐水6秒。如预期的,在肺动脉中的最大增强以后5-6秒,观察到升主动脉峰中的增强。同样地,升主动脉隔室中的增强水平比预期的肺动脉中的增强水平低约100HU,因为快速推注在肺动脉干和升主动脉之间被稀释。通过观察在50秒时肺动脉中的第二对比增强峰和接近60秒的第二峰在升主动脉中的出现,可以理解造影剂的全身再循环。
进行了另一个实验,其中进行了一系列模拟。对于一位模拟的患者,造影剂注射参数在集合中保持恒定,以证实所述模型的重复心输出量对对比增强的影响的能力。注射参数是:370mgI/ml浓度造影剂,在5ml/s注射10秒(体积=50ml),随后是在5ml/s的30ml盐水冲洗相。以1L/min的增量,将虚拟患者的心输出量从低值(3L/min)调至高值(8L/min)。该实验假定,可以独立于血液体积来操纵心输出量。就每个心输出量值而言,虚拟患者的血液体积保持恒定(使用与前一实验相同的参数:170磅,68英寸,35岁男性)。使用比在前一实施例中更短的注射持续时间来避免在低心输出量值时的任何再循环现象。
记录升主动脉中的峰值对比增强水平和最大对比增强的时间,并图形地显示在图4中。我们预期会发现,随着心输出量增加,造影剂快速推注在血管领域中的到达时间减少。同样地,随着心输出量增加,我们预期,血管结构中的对比增强会下降,正如在以前的研究中理论地和经验地证实的。
还使用虚拟的35岁男性(高68英寸,体重170磅),在MATLAB/Simulink中模拟了所述混合模型。在该组模拟中使用3个不同的注射方案:370mgI/ml造影剂,在5ml/s持续10秒(体积=50ml),随后在5ml/s使用盐水驱逐剂6秒(体积=30ml),盐水驱逐剂流率为2.5ml/s,持续6秒(体积=15ml),或没有盐水驱逐剂。得到的时间增强曲线显示在图5中。所述模拟证实,对于其中在造影剂快速推注的相同流率注射盐水冲洗的方案,峰值主动脉增强增加了14%。达到峰值增强的时间也延迟了0.8秒。对于其中在造影剂快速推注的流率的一半注射盐水的注射中,峰值增强的增加比不施用盐水冲洗时大7%。峰值增强时间偏移了0.6秒。Bae模型不能重复这些结果,因为Bae模型的公式不允许考虑在造影剂快速推注施用以后的盐水驱逐剂或稀释剂冲洗相的影响。所述结果证实了模型响应与其他人在动物模型和人类中证实的结果的定性一致。
如上所述,将混合PBPK模型结果与来自Bae模型的结果相对比。来自在Bae,K.T.,J.P.Heiken,和J.A.Brink,Aortic and hepatic contrastmedium enhancement at CT.Part I. Prediction with a computer model.Radiology,1998.207(3):第647-55页中公开的实验组群的集合数据促进了性能对比。基于造影剂注射方案,将患者登记成3组。在研究中使用了3个对比方案:二相低流率方案、单相低流率方案和单相高流率方案。在该研究的任何组中,没有使用盐水冲洗相。
使用从单相高流率组群得到的数据,评价新模型。在该组中,27位受试者的平均体重是177磅,值范围是44.1-135.0磅。对比方案包括:在5.0ml/s注射320mgI/ml造影剂,持续25秒(体积=125ml)。没有报告关于该组群的身高、性别或年龄组成的数据,但是在模拟他们的模型时,作者使用68英寸的平均身高和男性性别。每隔15秒对每位受试者在中间腹主动脉处进行单水平扫描120秒,然后每隔60秒扫描1次,直到300秒。在将1cm2ROI放在腹腔干处的腹主动脉上并将另一个放在肝实质上以后,由研究人员建立TEC。
在MATLAB/Simulink中执行混合PBPK模型。因为收集实验结果至300秒,造影剂具有足够的时间来在血管系统中再循环。将清除期加入所述模型中,以确保在时间增强曲线或TEC中的适当的下坡/再循环动力学。在Bae模型的模拟中设置典型的肾小球滤过率(GFR)(标称50-70ml/min),将清除期设定为模拟的肾动脉中的血液体积流率的19%。
[180]所述混合模型的许多实施例不含有肾隔室,在这样的实施例中模拟清除的机制可以是例如通过肺隔室。在对比公开的Bae实验数据的模拟中,使用0.08sec-1的血液-组织转换系数(kLUNG_BT)和0.1秒的隔室清除值(CLLUNG_T)来接近60ml/min的清除率,因为就该问题而言,没有有意义的生理学过程具有短于1秒的时间量程。太短的时间步会导致无效模拟,因为需要过多的时间来递送与用更长的时间步所得到的等效的结果。在0.01和0.001的时间步进行的模拟揭示了预测的增强没有差异。
使用来自Bae研究的集合患者数据的混合模型模拟的结果显示在图6中,紧靠来自集合的增强数据的结果。使用固定步长解算机(Runge-Kutta)和0.1秒的步长时间,进行所述模型的模拟。为了更好地显示与实验数据的对比,将增强数据降低采样。
在Bae,K.T.,J.P.Heiken,和J.A.Brink,Aortic and hepatic contrastmedium enhancement at CT.Part I.Prediction with a computer model.Radiology,1998.207(3):第647-55页中的高流率组的平均峰值增强被报告为313.7HU,而新模型产生了具有317.3HU的峰值增强的主动脉增强曲线,差异百分比为1.15%。Bae整体模型产生了具有321.3HU的峰值增强的增强曲线,差异百分比为2.42%。达到所述模型预测的峰值增强的时间是31秒,达到Bae整体模型预测的峰值增强的时间是31秒,实验数据具有32秒的平均峰值增强时间。
在Bae,K.T.,J.P.Heiken,和J.A.Brink,Aortic and hepatic contrastmedium enhancement at CT.Part I.Prediction with a computer model.Radiology,1998.207(3):第647-55页中呈现的验证结果没有评估所述模型在模拟单个对比增强测线中的实用性。为了确定Bae整体模型和本文的PBPK模型的预测单个对比增强测线的能力,使用来自CTA成像临床试验中的受试者的单个患者属性(身高、体重、年龄、性别)来呈现模拟结果,以将两种模型参数化。然后使用来自那些患者的成像数据的对应的增强测线(从临床图像提取的数据),对比2种模型的输出。
使用3种度量来评价2种模型的预测能力:均方根误差(RMSE),最大增强差异百分比(PDME),和在Bae,K.T.,J.P.Heiken,和J.A.Brink,Aortic and hepatic contrast medium enhancement at CT.Part I.Prediction witha computer model.Radiology,1998.207(3):第647-55页中描述为“增强差异指数(EDI)”的度量,所述EDI是,从临床成像数据集得到的模拟的和实验的增强曲线之间的差异的总和除以实验数据的曲线下面积:
其中下标“Sim”和“Emp”分别代表模拟的响应数据和实验的数据。
使用CT血管造影术研究数据集来对比所述混合模型和Bae整体模型的性能。该研究的目的是,使用在肺动脉干水平的定时快速推注扫描的特征,评价患者特异性的、优化的造影剂递送算法。在研究中登记了70位接受常规冠状动脉CTA(cCTA)的受试者,但是在本文呈现的对比实验中使用了来自总数中的20位的数据,因为尽可得到该数目的测试快速推注增强数据和降主动脉数据,这是由于不完全的轴向CT数据集贮存和转移。
在该分析中使用的20位受试者的人口统计数据总结在表3中。在该样品中的平均体重是美国成年人的平均值(参见Flegal,K.M.,等人,Prevalance and trends in obesity among US adults,1999-2000.JAMA,2002.288:第4页),且平均年龄是通常具有cCTA指征的人群的代表年龄。
表2.概括的人口统计数据平均数据,具有标准差
使用标准的扫描参数,用Dual Source CT扫描仪(Definition DS,Siemens Medical,和Malvern PA),扫描研究中的所有受试者。给所有受试者施用20ml测试快速推注(300mgI/ml造影剂),随后施用30ml盐水冲洗,流率为5ml/s。在开始造影剂注射以后大约5秒(由扫描仪软件引起的变动),开始定时快速推注扫描,并每隔2秒获取单水平扫描,直到升主动脉中的峰值增强以后大约5秒。扫描仪操作员慎重地停止数据获取,以减少不适当的辐射向受试者的暴露,没有得到获取数据的一致持续时间。
基于由扫描仪操作员在肺动脉干和升主动脉中画出的目标区或ROI,由扫描仪软件建立TEC。研究软件使用在这2个区域中达到峰值的时间和峰值增强来计算患者特异性的造影剂流率、体积、造影剂/盐水混合物和个体化的扫描延迟,用于诊断性CT扫描。总结的操作特异性的结果呈现在表3和表4中。许多因素会影响CTA检查的扫描持续时间,包括扫描仪设置、患者心率和对解剖学成像的持续时间。通过在公开的PCT申请号WO 2009/012023中描述的研究软件,计算在该数据集中的扫描延迟。
所述研究软件计算的对比方案由3个相组成:仅造影剂的部分,造影剂和盐水混合物相,和最后的仅盐水“冲洗”相。在表4中报告的相2的体积表示该相中的造影剂的体积,而相3体积是在冲洗中递送的盐水的量。因为在相2中稀释了造影剂,在表4中报告的浓度不是储备溶液的浓度(300mgI/ml)。
表3.来自临床数据集的相关扫描参数,前20位受试者。平均值给出了标准差。
表4.来自临床数据集的总结的对比方案统计学,前20位受试者。值是平均值±标准差。
使用半自动化的分段软件工具来从20位受试者成像集合中提取时间增强数据(螺旋成像数据来自诊断扫描)。每个轴向CT图像间隔0.375mm。获取数据块之间的时间在0.75和1秒之间变化。一个实施例数据集(受试者1)显示在图7中。该特定获取的扫描延迟是23秒,因而时间轴在23秒开始。因为增强数据是时间和空间的函数,将跨数据集的对应的空间尺度叠加在图的x-轴上。因为PK模型会产生关于在离散的“隔室”中的时间的预测,本文采用的对比是针对从20位受试者的主动脉提取的图像数据的头部侧位。
为了确定HU值和造影剂的血液-血浆浓度之间的关联,在临床研究之前,在使用的扫描仪上进行了使用稀释量的造影剂的实验。将稀释的造影剂放入不透射线的容器和扫描仪中,它们使用与在临床试验中相同的参数。通过线性回归,计算出HU和mgI/ml之间的换算因子为27.1HU/(mgI/ml)。该常数用于混合模型和Bae模型模拟中,以将碘的血浆浓度转换成HU。
为了使用临床成像数据提供所述混合模型相对于Bae模型的对比,在MATLAB/Simulink中实现了Bae整体PK模型。在构建所述模型时,使用在Bae,K.T.,J.P.Heiken,和J.A.Brink,Aortic and hepatic contrastmedium enhancement at CT.Part I.Prediction with a computer model.Radiology,1998.207(3):第647-55页中描述的、并在美国专利号5,583,902中进一步详细说明的模型的参数。固定时间步解算机(Runge-Kutta)每隔0.1秒执行所述模型。使用方程式27、28和29计算出单个受试者人口统计数据(身高、体重、性别和年龄),用于计算每个模拟的受试者的心输出量和总中枢血液体积估测值。如公开的,因为Bae模型不允许模型化盐水冲洗相,在模拟所述模型中使用的对比方案仅包括造影剂。使用27.1HU/(mgI/ml)的换算因子,将血液-血浆浓度转换成HU值。
使用在4.1.2部分中呈现的参数和方法,模拟混合模型。与Bae模型一样,受试者人口统计数据和在方程式27、28和29中定义的关联会提供心输出量和血液体积的值。在表1-3中的参数值用于构建所述模型。使用MATLAB/Simulink(R2008b)来模拟每位患者的模型,其中使用患者参数的临床数据和诊断相造影剂注射方案,包括盐水冲洗相。与Bae模型一样,使用27.1HU/(mgI/ml)的换算因子,将血浆浓度转换成HU值。
将每位受试者的预测HU值降低采样,以匹配来自成像数据的增强曲线的时间分辨率(1秒/样品)。在计算对比度量中,仅使用了跨成像数据的时间段的预测的增强曲线部分(MSE、PDME、EDI)。
所有20组人口统计、操作和成像数据足以允许模拟所述模型,并允许计算度量。在图8中,与来自成像数据的对比增强测线一起,绘制了使用受试者11的人口统计和操作数据的实施例模拟输出。在每个成像数据点上的误差棒是,通过分析在每个获取位置处在主动脉中的所有像素所测得的对比增强的标准差。在CT扫描仪、造影剂流动动力学、图像处理机制和结构活动中的固有成像噪声会促进在每个点处的噪声。
使用得自受试者11(60岁女性,体重280磅,身高68英寸)的数据的混合模型,产生了增强曲线,所述增强曲线与从该患者的主动脉成像数据获取的增强曲线较好地匹配(RMSE为15HU)。该患者的结果显示在图8中,因为预测的增强曲线与测得的数据非常好地匹配,也因为更重的患者呈现对常规造影剂给药方案的挑战。在4.1ml/s,递送了共61ml300mgI/ml造影剂,扫描持续时间为11秒。
Bae模型产生适当的增强曲线,但是会使实验数据反冲20-40HU。所述混合模型预测出比Bae模型更高的峰增强值。临床数据的线性外推提示,实际的(假定在体内测量主动脉的对比增强多10-20秒的能力)峰值增强值应当大于Bae模型的提示值,这与混合模型预测相一致。在图9A和9B中,给出了对比来自所述模型的增强预测的更多实施例和2组成像数据。图9A的对比是针对受试者6,她是37岁女性,体重110磅,身高61英寸,心率80bpm。递送77ml 300mgI/ml造影剂,流率为5.9ml/s,扫描持续9秒。图9B的受试者8是53岁男性,体重223磅,身高70英寸。他的心率平均值为58bpm。在4.1ml/s,注射74ml 300mgI/ml造影剂,他的扫描持续时间为12秒。
在某些情况下,Bae模型和混合模型预测的对比增强测线不同于实验数据。就受试者6而言,Bae和混合模型都预测比实验数据更小的增强最大值200HU。就受试者8而言,所述混合模型预测比实验数据的最大值小的最大增强79HU。Bae最大预测比实验数据的最大值小111HU。Bae和混合模型都预测比实验数据更小的增强最大值200HU。
对比混合模型和Bae模型的预测的统计结果呈现在下面的表和图中。就每位受试者而言,所述患者和操作数据用于两种模型中。在每种情况下,将所述模型的主动脉隔室增强输出与成像数据集相对比。在均方根误差和增强差异指数的计算中,包括从与扫描延迟对应的时间点开始并持续至扫描结束的模拟数据区段。结果列表显示在表5中。然而,所述混合模型模拟的3个对比度量更小,曼-惠特尼U检验没有揭示所有测量的中位值之间的显著差异(RMSE-U=455,p=0.229;PDME-U=445,p=0.351;EDI-U=470,p=0.064)。
表5.与临床数据相对比的模型总结结果
3个比较试验的箱线图揭示数据的偏态分布和中位值的等价。用延伸至比中位值大或小1.5倍四分位距(IQR)的数据点的短线,绘制在图10中的箱线图。在图10中的十字符号指示偏离中位值超过1.5*IQR的数据点的轮廓。
所述数据的另一种可视化作为在图11中的一系列散布图给出。图表的目检揭示,受试者6、8、12、16和19的结果显示出模拟性能的最大变化。与Bae数据相比,受试者15-20的混合模型预测值显示出与临床数据稍微更好的一致性。所述混合模型在20位受试者中的14位中产生更低的RMSE值,在20位受试者中的13位中产生更低的PDME结果,在20位受试者中的17位中产生更低的EDI。
与Bae等人的公开的整体的生理学基础上的药代动力学相比,本文中的混合模型具有降低的阶,是可离散化的,会明确地模型化造影剂传输延迟,会模型化在CTA操作过程中常规地进行的盐水冲洗注射的效应,并模型化在外周静脉循环中的随时间变化的血液-血浆/造影剂相互作用效应。如上所述,所述混合模型的输出与Bae等人公开的临床数据的对比,会显示有利的一致性。所述混合模型与从螺旋CTA数据导出的对比增强测线的进一步对比,揭示了与Bae模型的紧密一致性,且在许多情况下,揭示拟合临床数据的更好能力。
下面描述了所述混合模型用于测试方法的用途,所述方法鉴别来自测试快速推注数据集的模型。在许多研究中,采用简化策略,而不是尝试拟合混合模型的所有参数。例如,可以进行该模型简化,以减少计算负担和数字预测的可靠性。
数据驱动的方法
在许多实施例中,本文所述的方法、系统和模型用于开发患者特异性的对比增强预测,其中使用从测试快速推注导出的对比增强数据。评价了2种技术,一种假定模型结构(参数的),另一种使用非参数的“黑盒子”方案。在所述基于模型的方案中,将模型简化策略应用于上述的混合药代动力学模型,以减少在鉴别模型参数时的计算负担,并克服随时间而变化的系统的建模挑战。
其它方案,包括参数的和非参数的方案,参见,例如,公开的美国专利申请号2007/0255135,公开的美国专利申请号2007/0282263,公开的美国专利申请号2008/0097197,公开的美国专利申请号20080097197,公开的美国专利申请号2010/0113887,公开的美国专利申请号2010/0030073,和公开的PCT国际专利申请号WO 2009/012023,它们的公开内容通过引用并入本文。
参数的(基于模型的)鉴别
在许多用于构建数据依赖性的、患者特异性的造影剂传输模型(具有上述的降阶形式的混合PBPK模型)的实施例中,在获取例如2条由注射小“鉴别”快速推注造影剂(10-20ml)引起的CT时间-增强曲线(在肺动脉干水平的相继低水平CT扫描)以后,估测模型参数。
将造影剂输入外周循环系统,通常通过在前臂静脉中或在肘窝中的血管内入口。几秒以后,造影剂快速推注到达右心。接着,右心室推动造影剂快速推注穿过肺循环。在此时,造影剂快速推注的运输由中枢循环参数(即心输出量)支配。在6-20秒以后,造影剂到达心脏左侧,并被左心室射入主动脉和冠状动脉脉管系统。通过将轴向CT获取定位在肺动脉干水平并每隔n秒在该水平获取CT图像,产生造影剂的运输的数字和图形描绘。
图12呈现了从经历CT血管造影术的患者获取的横截面图像数据的轴向“切片”。在外周注射造影剂的小剂量测试快速推注(20ml,350mgI/ml),并每隔2秒在肺动脉干水平进行扫描。随着造影剂流过解剖学,血管结构中的图像增强或量度增加。在获取轴向图像以后,扫描仪或离线处理软件在每个图像上的目标区(ROI′s)内进行表面积分,并相对于时间绘制平均衰减值。标准的临床实践是,在造影剂注射开始以后5-10秒开始获取,并在标准管电压(120kV)、但是在低管电流-时间乘积(10-30mAs)获取图像20-30秒。
人受试者的典型时间增强曲线或TEC显示在图13中。将目标区(ROI)标志物放在受试者的肺动脉干和升主动脉的上面,如图12所示。
没有尝试使用从测试快速推注TEC可得到的有限持续时间数据来拟合本文中的混合模型的所有参数。相反,采用模型简化策略,其中在造影剂的首次通过过程中,考虑所述混合模型的最相关的隔室:外周静脉室、右心隔室、肺隔室和左心/主动脉隔室。下述共识证实了该方案:在造影剂的首次通过过程中进行CTA成像,且通常在造影剂注射开始以后30秒时结束,这归因于扫描仪的短扫描获取时间。当在系统中存在反馈时,也存在对统计鉴别技术的固有挑战。造影剂再循环被视作所述混合模型中的一类反馈。
图14是用于参数估测的混合模型子系统的一种图解。将PER和RH隔室收缩成1个传输函数He1,其中造影剂注射输注是输入函数uinj(n),且测量的TEC(图13中的ROI1)是用于参数估测的混合模型子系统的一种图解。将PER和RH隔室收缩成1个传输函数He1,其中造影剂注射输注是输入函数uinj(n),且在肺动脉干中测量的TEC(图13中的ROI1),yRH(n),是输出信号。使用取决于扫描仪构型的缩放系数,将yRH(n)转换成血浆浓度,并通过校准进行测定。因为在所述模型中假定线性,描述从注射部位至肺动脉的造影剂快速推注的传输延迟的运输延迟可以定位在He1块内的任意位置。同样地,用传输函数He2,将LUNG_B和LH子系统组合进一个子系统。He2块的输入函数是yRH(n),输出信号是从ROI(图13中的ROI2)导出的测量的浓度TEC yLH(n),所述ROI在肺动脉干的测试快速推注扫描过程中放在主动脉中,yRH(n)是输出信号。使用取决于扫描仪构型的缩放系数,将yRH(n)转换成血浆浓度,并通过校准进行测定。因为在所述模型中假定线性,描述从注射部位至肺动脉的造影剂快速推注的传输延迟的运输延迟可以定位在He1块内的任意位置。同样地,用传输函数He2,将LUNG_B和LH子系统组合进一个子系统。He2块的输入函数是yRH(n),输出信号是从ROI(图13中的ROI2)导出的测量的浓度TEC yLH(n),所述ROI在测试快速推注扫描过程中放在主动脉中。
通过在2种不同的操作中进行参数估测,减少计算复杂性。除了减少计算复杂性以外,当使用输入/输出数据鉴别更少的参数时,可以实现参数差异的理论减少。此外,使用一对输入/输出数据鉴别出的参数总数的减少,有助于防止过度拟合(overfitting)。
在初步的参数估测过程中在HE1系统中仅包括RH和PER隔室,会导致实验数据的不收敛和差拟合。通过在外周静脉隔室和右心隔室之间引入“中间”隔室来增加He1系统的阶数,会导致参数估计值的收敛。该中间隔室也会提供模型化周边隔室内的注射流率效应的方式,因为在周边隔室之间的质量通量是由新的体积流率术语QPER驱动,该术语QPER既不是注射流率,也不是进入右心的血液的流率。所述模型结构呈现在图15中,它显示了:uinj(t),即向外周静脉中的造影剂快速推注,和yRH(t),即在肺动脉中测量的TEC(并从HU转换成浓度单位)。
He1的连续时间状态空间公式是:
BHe1=[1 0 0]T
DHe1=[0] (31)
yHe1(t)=CHe1xHe1(t) (32)
并且,当测量的数据可从在右心结构(诸如肺动脉)中收集的CT数据得到时,调用yHe1(t)=yRH(t)。穿过外周静脉循环的造影剂的运输延迟是TPER。用施用流率(Qinj)、测试快速推注浓度(Cinj)和造影剂快速推注持续时间TDUR,将方程式32的输入函数参数化:
uinj(t)=Qinj(t)Cinj[u(t)-u(t-TDUR)] (33)
其中u(t)是单位阶跃函数。
在方程式31的公式中省略造影剂注射流率(Qinj),因为假定,当在测试快速推注造影剂体积以后进行盐水快速推注冲洗时,获取识别数据集。在测试快速推注以后的盐水输注会维持周边隔室中的流率(QPER)10-14秒的持续时间(造影剂和盐水的体积除以流率)。典型的诊断快速推注具有10-14秒的持续时间(造影剂),所以假定是,在诊断注射过程中的外周流率与在测试快速推注过程中相同。
因为在不同的时间步通过单水平扫描获取CT数据,在识别步骤中使用的模型被离散化。通过线性系统理论的标准转化,得到在方程式31中定义的连续状态空间系统的离散化:
下标“D”表示离散化的状态空间矩阵形式,Δ是由单水平扫描之间的时间限定的取样间隔。该方案会产生下述状态空间矩阵:
DHe1D=[0] 35
所述模型的第二个子系统是肺和左心隔室的组合,在图14中表示为He2。因为首次通过循环效应在预测CT血管造影术的对比增强时占优势,可以忽略造影剂从血管区向肺组织的运输。正如在所述混合模型的衍生化中,穿过心肺回路的造影剂的传输延迟汇集成运输延迟TLUNG。分离的He2子系统的图形描绘显示在图16中。
He2的连续时间状态空间动力学是:
BHe2=[QCO 0]T
DHe2=[0] (36)
yHe2(t)=CHe2xHe2(t) (37)
并且,当TEC数据可从左心结构(诸如升主动脉)得到时,yHe2(t)=yLH(t)。将方程式34应用于系统方程式,会产生A和B矩阵的离散时间表达(C和D矩阵等同于连续时间形式),其中Δ表示取样时间:
DHe2D=[0] (38)
2个离散系统充当本文开发的数据驱动的参数估测技术的实施例的基础。向第一个离散系统(He1D)中的输入是取样的注射输入信号,即uinj(n)(方程式33的离散形式),输出是在肺动脉干中测得的TEC,即yRH(n)(图13中的ROI1)。向第二个离散系统He2D中的输入是yRH(n),输出向量是从升主动脉测得的TEC,即yLH(n)(图13中的ROI2)。图17显示了2个离散系统的块状图和所述系统的各个输入和输出。在参数估测中使用的2个离散的状态空间系统是:
要在第一个系统He1D中估测的参数是:VPER1、VPER2、VRH、QPER、QCO和运输延迟TPER。要在第二个系统中估测的参数是:VLH、VLUNG、QCO和穿过肺系统TLUNG的运输延迟。
最大似然估计值
在许多从在将小剂量造影剂快速推注或物质注射进患者中以后获取的扫描数据估测参数的实施例中,参数估测发生在2个不同的步骤中。通过使用2个测量信号yRH(n)和yLH(n),并系列地估测2个系统,会减少收敛时间、计算负担和参数变异性,因为通过数字优化,在一个单独的识别步骤中需要识别更少的参数。例如,不是使用取样的测试快速推注增强数据估测10个参数,而是使用已知的输入函数uinj(n)和输出信号yRH(n)(其可能仅具有10个数据点)识别He1D的6个参数,然后在另一个单独的参数估测步骤中识别He2D的仅4个参数(使用yRH(n)和yLH(n)作为输入/输出数据)。
首先,使用得自测试快速推注的扫描数据和在肺动脉和升主动脉中测量的得到的TEC(yRH(n)和yLH(n)),将“He2”系统参数化。接着,使用信号uinj(n)和yRH(n)、输入注射函数和在肺动脉中测量的TEC,将“He1”系统参数化。在估测参数以后,使用方程式39和40,构建该患者的单个模型。通过将升主动脉中的预测的增强(使用测试快速推注)与测量的数据相对比,确定拟合优度。
因为大量噪声源会影响测量过程(包括零阶取样过程),可以假定,例如,误差在系统中正态分布。在高斯分布误差和独立噪声过程的假设下,最大似然估计(MLE)可以用于估测未知的模型参数。本领域技术人员显而易见,其它方法可以例如用于估测未知参数。当估测He1D和He2D的参数时,在MLE中使用的代价函数是测量的测试快速推注TEC数据y
RH(n)和y
LH(n)与那些结构中估测的测试快速推注响应
和
之间的平方化差异的总和:
其中参数估测向量(“hat”符号,^,表示估测的参数):
公认的是,当噪声过程是高斯(Gaussian)时,最大似然估计等同于最小二乘法拟合。在方程式42和42中的方差术语
没有输入代价函数,且可以在事后估测,基于得自测量的TEC的噪声估测。在方程式41和方程式42的总和中,N是可从测试快速推注增强数据得到的TEC样品的数目,指数变量i的范围跨离散时间样品。
在矩阵符号中,将问题阐述为,通过估测的参数的函数,使残差的平方和最小化:
]
其中
是在给出p-值参数估计
和
(p=6和p=4)和输入的情况下估测LH和RH信号的n-向量。向量
和
是He1D和He2D系统的估测的残差:
其中η(n)表示模型和测量不确定度。
从残差向量相对于参数向量的一阶导数,构建出nxp雅可比矩阵
和
应当指出,粗体的
符号用于雅可比矩阵,而V(非粗体)用于体积。雅可比矩阵的第i个和第j个元素是:
参数向量的最大似然估计呈正态分布,协方差为∑=σ2(VTV)-1,其中方差的估测器是[55]:
参数估计的方差还可以表示为费希尔信息矩阵的倒数(黑塞矩阵的倒数),其经常可从数字计算包(诸如MATLAB的优化工具箱)得到。
参数偏倚(它是估测的和真实的参数之间的差异)的近似表达,可以以雅可比矩阵的方式来表示:
因为由方程式39和方程式40定义的模型的参数是非线性的,代价函数的最小化是非线性的最小二乘法(NLS)问题。存在众多用于求解NLS问题的迭代技术。该问题的一个优化方案可以提供对参数的非负性约束。并且,对参数估计施加生理学约束是合理的。因为该问题是受约束的最小化,一维搜索技术(如Gaus-Newton)、单一的或混合的技术(如经典的Levenberg-Marquardt)不可适用。为此原因,使用子空间置信域方法(其基于内反射Newton方法,如在MATLAB优化工具箱(nonlinsq)中实现的)来估测方程式43中的参数向量。该算法使用预条件共轭梯度(PCG)技术计算解算机的每个迭代处的近似解。在求解过程中,使用有限差分逼近来构建雅可比矩阵和黑塞矩阵。
在下述的实验中,运行最小化,在代价函数评价上的收敛判据为10
-5,在参数上的耐受性为10
-5。设置的函数评价的最大数目是400,允许的在每个求解步骤处的迭代的最大数目是500。在实验中使用的
和
的参数边界呈现在表6中。
表6.对估测的参数的上和下边界或约束
在开始估测之前,将参数初始化在上和下参数边界的50%中点处。为了研究解算机对初始条件选择的灵敏度,每个参数估测执行5个解算机运行。选择具有最佳拟合输出曲线的初始化向量作为起始估测。当执行对比增强的正向预测时,使用He1D系统的估测的输出向量
作为He2D系统的输入。
参数估测器评价方法。
为了证实参数估测技术的有效性,进行了3个不同的实验。前2个实验使用合成的数据来确定估测器的性质(偏倚、方差)和它对截短的输入/输出数据及数据的暂时区域的变动的稳健度。然后使用回顾临床数据来确定估测器在预测人受试者的主动脉对比增强中的性能,这基于在肺动脉干和升主动脉中测量的定时快速推注扫描TEC。
模型至模型对比。为了表征估测器,使用由方程式39和40(由估测器使用的相同方程式)定义的系统,建立合成的数据集。因此,所述混合模型没有产生在这些实验中使用的数据。该方法的目的是,分离独立于测量数据的估测器的性能。这类对比在文献中已经被称作“模型至模型”对比,并用于研究独立于其它因素(诸如模型精确度,采样噪声和噪声过程中的变异性)的估测器的性能。所述模型至模型术语表示,用于估测参数产生的数据的相同模型。它也可用于直观化解空间和评估估测器的性能,因为已知“真”参数
(就i=1、2而言)。使用在表7中指出的参数值,计算数据集。
表7.在模型至模型模拟和对比中使用的参数值
将方程式33的输入函数uinj(t)参数化,造影剂浓度为370mgI/ml(Cinj),持续时间为5秒(TDUR),流率为5ml/s(Qinj)。使用0.10秒/样品的时间步,在MATLAB中进行模拟。在1秒/样品、2秒/样品和5秒/样品的速率,将模拟的TEC降低采样,以模仿CT扫描仪对TEC的采用作用。
还用“真”参数进行模拟,所述“真”参数用合成的测量噪声(AWGN)校正过,所述AWGN具有0、1和2HU的标准差(在方程式46和47中的η(n)),以产生偏倚和方差的估测。对于噪声水平和样品时间的每种组合,进行30次参数估测。对于每次模拟运行,将初始参数估计改变标称值的25%。计算5个参数估计,并通过选择产生在得自参数估测的预测增强曲线和由模型
产生的增强曲线之间的最小剩余均方差的集合,来确定最佳的参数集合。
最后,使用真参数值进行无噪声模拟,所述真参数值以2%增量变化±20%,以产生代价函数JRH(θHe1D)和JLH(θHe2D)的图形描绘。这些模拟也允许分析估测器对初始条件选择的灵敏度。执行5次模拟,其中使用在表7中的值的25%内随机变化的初始条件。记录并分析每个初始条件集合的预测的测试快速推注和模拟的测试快速推注增强之间的均方根误差。
使用混合模型数据的估测器性能。使用混合模型数据,在有噪声和对截短的测试快速推注增强数据的灵敏度存在下,测定估测器的性能。使用上述的参数和方法来模拟混合模型,其中模拟的受试者人口统计数据和在方程式27、28和29中定义的关联会提供心输出量和血液体积的值。在表1-3中的参数值构建所述模型。使用MATLAB/Simulink(R2008b)模拟每个患者的模型,其中使用患者参数的临床数据和诊断相造影剂注射方案,包括盐水冲洗相。使用相同的参数集合来从方程式9和11建立合成的、离散的测试快速推注数据
和
其中上标(Test,H)表示从所述混合模型建立的模拟的测试快速推注响应。用于上述模拟的患者人口统计数据会参数化在方程式15中阐述的模型系统。将固定的流率和造影剂体积注射信号(uDiag(n))递送给每个模拟的“受试者”,以产生模拟的增强曲线
使用下述值,将注射方案u
Diag(n)参数化:Qinj=5ml/s,Cinj=350mgI/ml,持续25秒,继之以盐水相,流率为5ml/s,持续时间为8秒。
使用高斯噪声向量进行参数估测,将0.1、0.25、0.5、1.0、2.0、5.0、10和20HU的标准差加入混合模型测试快速推注数据中。使用在4.0RMSE中产生的度量、最大增强差异百分比(PDME)和增强差异指数(EDI),评估预测的增强相对于通过20个数据集的混合-模型模拟的“真”增强的性能。
对技术人员而言,进行研究来在峰值对比增强以后不久停止从测试快速推注扫描收集数据,是常规实践。停止数据获取的原因包括:担心向患者过度暴露辐射,和希望进行短时间研究。因为测试快速推注数据记录可以含有在峰值增强以后的仅少量数据样品,患者特异性的造影剂递送算法必须在有截短的数据集存在下起作用。
使用以模拟的混合模型建立的截短的测试快速推注数据向量进行实验,以测定截短的测试快速推注扫描数据对MLE参数估测器的影响。在20、25、30和35秒的持续时间,且不存在可加噪声,截短向量。接着,在25和35秒,并加入噪声,截短模拟的测试快速推注TEC。使用0.1、0.25、0.5、1.0、2.0、5.0、10和20HU的标准差,建立AWGN。在所有情况下,
和
的第一个数据点是开始注射以后5秒。选择25和35的截短时间作为具体的试验点,因为在临床实践中可以预期25秒获取持续时间,且35秒会确保捕获造影剂峰和在所述峰以后的许多样品(就大多数操作而言,在左心隔室中)。
使用回顾地收集的临床数据的估测器性能。在南卡罗来纳医科大学(Somatom Definition DS,Siemens Healthcare Malvern PA),针对在IRB批准的临床试验过程中收集的临床数据,评估估测器性能。数据集和用于收集数据的方法如上所述。使用在临床试验过程中收集的测试快速推注TEC来导出参数估计,并使用半自动化的主动脉对比增强数据来对比估测器的预测输出,与上述的方法相同。
在临床试验过程中,扫描仪操作员通过将ROI放在肺动脉干和升主动脉上来获取测试快速推注扫描数据。在造影剂注射开始以后5秒,开始数据收集。在观察到升主动脉中的对比增强峰以后2-4秒,由操作员停止单水平扫描获取。扫描仪软件然后处理数据,并产生TEC,所述TEC输出为数据文件,并保存。如上所述提取诊断扫描增强测线数据。
实验方法的概要是:(1)提取测试快速推注TEC数据;(2)使用MLE技术,进行参数估测;(3)使用得自临床数据集的诊断注射方案和得自(2)的鉴别参数,产生预测的对比增强;(4)从临床数据集提取增强曲线;(5)将预测的对比增强与临床数据(从扫描延迟至扫描获取结束)进行对比。
将得自数据驱动的估测器的预测输出与使用RMSE、PDME和EDI的实际的临床数据相对比。
使用MATLAB(r2008b)、优化工具箱(v4.1)和Simulink(v7.2),进行所有模拟和分析。
使用模型至模型对比的估测器性能。得自所述模型-模型模拟的结果显示在表9至表12中。在表9中,列出了心肺系统He2D的每个参数的得到的平均估计值偏倚。用于衡量估测器的性能的拟合优度判据是预测的
和估测的yLH之间的均方差。除了每个参数的偏倚百分比以外,还将跨30次模拟的MSE均值(和MSE的平方根)制表。在表10中,总结了He1D子系统参数估测的偏倚结果。
在He2D子系统中,参数估测偏倚独立于可加噪声而增加,在0.5和1秒/样品的采样周期时的Qco和TLUNG参数除外。所述拟合的MSE似乎独立于采样周期,因为它作为所有3个采样周期的噪声σ的函数类似地增加。当数据采样周期超过1秒/样品时,He1D子系统中的参数估测偏倚增加。在估测器的预测值和真值之间的残差的MSE随着噪声贡献增加而增加,采样周期大于1秒/样品。
表9.θHe2D的参数偏倚的总结。值是得自标称值的偏倚百分比。MSE是估测的和实际的增强信号yLH的残差的均方差。
Ts |
σHU |
VLH |
VLUNG |
QCO |
TLUNG |
MSE |
0.5 |
0 |
-2.2 |
9.6 |
3.7 |
-8.1 |
0.0 |
|
1 |
-6.6 |
11.6 |
2.5 |
-7.2 |
1.0 |
|
2 |
-10.5 |
18.1 |
4.3 |
-5.7 |
4.2 |
1 |
0 |
-5.1 |
8.2 |
1.7 |
-15.9 |
0.0 |
|
1 |
-4.5 |
14.3 |
5.2 |
-15.4 |
1.0 |
|
2 |
-9.2 |
23.0 |
7.0 |
-14.3 |
4.1 |
2 |
0 |
-14.9 |
18.4 |
2.9 |
-29.8 |
0.0 |
|
1 |
-13.7 |
16.7 |
2.8 |
-29.2 |
1.0 |
|
2 |
-17.3 |
32.2 |
9.9 |
-25.6 |
3.7 |
表10.θHe1D的参数偏倚的总结。值是得自标称值的偏倚百分比。MSE是估测的和实际的增强信号yRH的残差的均方差。
Ts |
σHU |
VRH |
VPER1 |
VPER2 |
TPER |
QPER |
QCOr |
MSE |
0.5 |
0 |
0.0 |
3.0 |
12.5 |
-0.4 |
0.4 |
10.0 |
0.0 |
|
1 |
-10.2 |
2.1 |
22.8 |
-0.4 |
0.3 |
9.9 |
1.0 |
|
2 |
-4.4 |
3.1 |
19.1 |
-0.4 |
0.4 |
10.1 |
3.8 |
1 |
0 |
0.0 |
2.7 |
12.5 |
-0.9 |
0.4 |
10.0 |
0.0 |
|
1 |
-4.5 |
2.3 |
17.4 |
-0.9 |
0.3 |
10.0 |
1.0 |
|
2 |
-4.5 |
1.3 |
16.2 |
-0.9 |
0.2 |
10.1 |
3.7 |
2 |
0 |
35.1 |
1.9 |
-90.0 |
-1.0 |
0.3 |
10.7 |
2.4 |
|
1 |
34.9 |
3.3 |
-90.0 |
-1.0 |
0.5 |
10.6 |
3.4 |
|
2 |
38.1 |
3.3 |
-90.0 |
-1.0 |
0.5 |
10.8 |
5.3 |
2个系统的参数估计的方差显示在表8和表9中。当噪声贡献的标准差是2HU时,对于所有取样时间,参数估测方差最大。运输延迟参数TLUNG的方差在所有参数中最大。
表8.He2D的参数估测方差的总结
Ts |
σHU |
VLH |
VLUNG |
QCO |
TLUNG |
0.5 |
0 |
4.14E-07 |
4.01E-07 |
2.02E-05 |
4.92E-03 |
|
1 |
2.01E-02 |
1.88E-02 |
9.50E-01 |
2.30E+02 |
|
2 |
3.38E-04 |
3.08E-04 |
1.55E-02 |
3.79E+00 |
1 |
0 |
4.20E-06 |
4.01E-06 |
2.02E-04 |
4.69E-02 |
|
1 |
9.80E-03 |
9.02E-03 |
4.56E-01 |
1.14E+02 |
|
2 |
1.52E-04 |
1.35E-04 |
6.92E-03 |
1.84E+00 |
2 |
0 |
4.87E-05 |
4.23E-05 |
2.16E-03 |
5.16E-01 |
|
1 |
5.25E-03 |
4.67E-03 |
2.34E-01 |
5.61E+01 |
|
2 |
7.00E-05 |
5.80E-05 |
2.88E-03 |
8.32E-01 |
表9.He1D系统的参数估测方差的总结
Ts |
σHU |
VRH |
VPER1 |
VPER2 |
TPER |
QPER |
QCOr |
0.5 |
0 |
6.23E-12 |
1.06E-09 |
7.83E-12 |
7.58E-08 |
5.96E-08 |
3.56E-10 |
|
1 |
3.83E-03 |
6.26E-01 |
4.55E-03 |
4.51E+01 |
3.55E+01 |
2.13E-01 |
|
2 |
5.85E-05 |
8.95E-03 |
7.09E-05 |
7.00E-01 |
4.99E-01 |
3.29E-03 |
1 |
0 |
3.50E-11 |
5.50E-09 |
4.39E-11 |
4.25E-07 |
3.09E-07 |
2.00E-09 |
|
1 |
1.88E-03 |
3.04E-01 |
2.32E-03 |
2.26E+01 |
1.70E+01 |
1.06E-01 |
|
2 |
2.83E-05 |
4.87E-03 |
3.47E-05 |
3.39E-01 |
2.71E-01 |
1.57E-03 |
2 |
0 |
2.22E-03 |
3.94E-01 |
3.59E-03 |
2.97E+01 |
2.19E+01 |
1.21E-01 |
|
1 |
3.14E-03 |
5.16E-01 |
5.08E-03 |
4.20E+01 |
2.87E+01 |
1.71E-01 |
|
2 |
1.87E-05 |
3.08E-03 |
3.05E-05 |
2.52E-01 |
1.70E-01 |
1.02E-03 |
为多种参数建立左心和右心估测器(方程式44和45)的无噪声代价函数的等值线图,并显示在图18和图19中。在每个图中,2个参数在标称值附近变化,而所述模型中的其它参数是固定的。得到的图是n-空间超平面的二维投影图,其限定了在指定范围内的参数空间。当延伸在解空间的最小值附近的椭球面时,这会指示,数字解算机具有会聚真实最小值和增加的方差的问题。所述等值线图会显示单个最小值。许多参数对会表现出在最小值附近的长的且狭窄的椭球面——尤其是He2D子系统的VLH、VLUNG对。
混合模型参数估测结果。下面显示了使用混合模型作为20位受试者的增强数据源时的参数估测实验结果。在所有20位受试者中,使用20ml的输入注射,在5ml/s注射350mgI/ml浓度造影剂,随后注射40ml盐水,在所述混合模型的左和右心隔室中产生测试快速推注增强数据。将AWGN(标准差=2.5HU)加入模拟响应向量中,在35秒截短测试快速推注增强信号。
[287]模拟的数据集的诊断注射是75ml快速推注(350mgI/ml,在5ml/s),继之以40ml盐水(在5ml/s)。如图20所示,在估测的和模拟的数据(每位受试者的扫描延迟和扫描持续时间来自临床数据集)之间的对比中,仅使用在扫描延迟和扫描结束之间的数据。对于所有受试者,在如上定义的RMSE、PDME和EDI的计算中,包括在2条垂直线之间的数据点。
在模拟的和估测的响应之间的平均均方根误差是7.78±4.40HU。20位受试者的2条曲线的最大增强差异百分比(PDME)的最大平均值是1.29±1.12%,平均增强差异指数是1.57±1.16%。
测试快速推注数据长度和可加噪声对估测器的性能的影响显示在图22中。对于所有受试者,使用相同的诊断注射(75ml 350mgI/ml浓度造影剂,在5ml/s,继之以40ml盐水),进行临床集合中的20个受试者的混合模型模拟。
在产生图22的实验中,可加噪声为零。随着测试快速推注数据向量下降,RMSE和PDME增加。意外的是,因为增加的患者辐射暴露,在峰值对比浑浊化出现以后许多秒,扫描仪操作员会获取测试快速推注数据。得自测试快速推注的升主动脉中的典型峰时间是17-24秒。在模拟组群中达到峰值的平均时间是21.1±2.1秒。当在20秒截短测试快速推注曲线时,在大多数情况下失去峰值增强。尽管如此,所述估测器仍然能够产生系统动力学的估测,具有12HU RMS误差。
图23阐明了当将0.1和20HU之间的可加噪声加入在25和35秒截短的测试快速推注增强数据中时,估测器的性能。就RMSE而言,存在大约3HU的一致分离,对于更长的测试快速推注增强曲线,存在更少的RMSE。最大对比增强的预测对测试快速推注测量数据向量的长度不太敏感。在临床测试快速推注增强数据上的典型噪声是在2-10HU范围内。25和35秒测试快速推注向量之间的PDME的差异低于0.5%。2种对比度量保持恒定,直到可加噪声σ超过2HU。然后,误差线性地增加。
使用临床数据的最大似然估计结果。本部分呈现了估测器的使用回顾临床数据集来预测对比增强的能力。对于每位受试者,将在肺动脉干和升主动脉中测得的测试快速推注TEC数据用作向MLE算法中的输入。在估测的模型(基于MLE拟合)和主动脉的诊断扫描对比增强之间进行对比。在进行对比时,使用得自临床数据集的注射方案和扫描参数,且在RMSE、PDME和EDI的计算中仅包括在扫描延迟和扫描结束之间的数据点(如图20所示)。
图24呈现了使用测试快速推注数据和MLE技术估测的所有20位受试者的结果。水平的虚线指示每个性能度量的样品平均值。受试者7、9和15在所有3类之间具有最大误差。在表13中显示了使用MLE的分析的数字结果,以及当使用混合模型(仅用受试者人口统计数据参数化)时的预测结果。MLE方法产生更低的平均RMSE、PDME和EDI(PDME是显著不同的,p<0.05)、所有3种度量的更小数据范围、和更低的值标准差。混合模拟的最大RMSE是142.3HU,而使用MLE的RMSE是70.5HU。
表13.对比MLE和混合模型预测结果的总结结果
图25描绘了2位受试者(6和8)的预测的对比增强和诊断对比增强。使用混合模型对这些患者预测的增强显示在图9中。它们的预测性能在所有20位受试者中最差。使用得自临床数据的测试快速推注TEC和MLE方法,实现了这些2位受试者的更好的对比增强预测(对于受试者6,39HU相对于160HU的RMSE;对于受试者8,21HU相对于57HU——MLE和混合结果)。
图26呈现的箱线图图示地对比了在20位受试者之间得自所述混合模型和MLE的结果。在每个箱中的水平线指示中位值数据值,垂直箱的边缘指示数据的第1和第3四分位数。
非参数识别
下面阐述了使用测试快速推注范例来预测对比增强的非参数识别技术实施例的方法和结果(模拟和临床数据)。引入的截短的奇异值分解方法,呈现了它对模拟和回顾临床数据的性能,并呈现了与MLE方法的对比。
如上所述,使用测试快速推注增强数据来参数化降阶药代动力学模型,开发了用于估测和预测在个体基础上的对比增强的参数识别方法。或者,可以使用非参数的(或不依赖于模型的)方案,其也依赖于得自测试快速推注和扫描的数据,但是,例如,需要仅得自一个扫描位置的数据。在几项研究中,使用升主动脉作为扫描位置。非参数方案的一个优点是,不需要假设,诸如模型结构和阶数。但是,例如可以做出假设,根本的系统动力学可以模型化为LTI系统。
所述非参数方法提出对比增强问题作为反对问题,其中已知输出(在施用测试快速推注以后得自升主动脉的图像数据)和输入(测试快速推注碘施用测线)。一种正则化方法,如下所述的截短的奇异值分解(tSVD),用于估测药物和心肺系统的脉冲响应(或残差函数)。本领域技术人员显而易见,其它方法可以用于估测脉冲响应。
可以为tSVD选择适合测量噪声且会平衡残差和解误差的截短指数,以生成希望的结果。如在本文的附录1(在说明书末尾处)中所讨论的,与Fleischmann/Hittmair傅里叶解卷积方法相比,适当的截短指数和tSVD对数据变动是更稳健的,因为在该方案中,仅使用具有固定截止频率的低通滤波器。
当建立预测性的对比增强技术时,用实施例证实了正则化技术(如tSVD)的益处。使用标准的线性最小二乘法(MATLAB的矩阵左除算子),产生在临床数据集(测试快速推注TEC数据)中的受试者7的估测的脉冲响应hsys。通过脉冲响应和从附录1的输入函数A1形成的下三角形Toeplitz矩阵的矩阵乘积,定义增强动力学:yLH=Uinjhsys。图27描绘了hsys的线性最小二乘解,其中使用从在5ml/s注射20秒的测试快速推注(350mgI/ml浓度造影剂)产生的左心的临床增强数据。显而易见由测量噪声导入的振荡。线性最小二乘法方法不能过滤噪声,这阻止了系统脉冲响应的稳健估测。
在图28中描述了使用单一最小二乘法的脉冲响应估测的对比,它是使用在本部分中开发的tSVD方法对受试者7的估测脉冲响应。使用tSVD建立对比增强预测的详细算法显示在下面的表14中。
表10非参数估测算法
研究并实现了在奇异值分解中用于选择截短指数k的2个方案:Koh等人的线性分段拟合方法,和Hansen等人的适当修正L-曲线标准。Koh TS,W.X.,Cheong LH,Lim CCT,Assesment of Perfusion by DynamicContrast-Enhanced Imaging Using a Deconvolution Approach Based onRegression and Singular Value Decomposition.IEEE Trans Med Imaging,2004.23(12):第1532-1542页.Hansen,P.C.,The truncated SVD as a method forregularization.BIT,1987.27:第534-55页。使用在正则化工具箱(在因特网上自由散布,tsvd.m)中提供的截短的奇异值分解方法,建立脉冲响应估测。Hansen,P.C.,Regularization tools:a MATLAB package for analysis andsolution of discrete ill-posed problems.Numer.Algorithms,1994.6:第35页。
如附录1所示,Koh方法通过用分段的线性函数逼近Picard图来计算截短指数。Picard图是傅里叶系数相对于指数变量的对数图,所述指数变量的范围是观察向量
中的数据样品的长度。该问题的线性方程式系统是:
它类似于线性系统的Ax=b的标准符号。傅里叶系数是左侧单个向量(得自奇异值分解的左侧单个矩阵的列)和观察向量的乘积的绝对值|ui TyTest|。通过下述求解,确定用于逼近Picard图的系数:
其中Y是傅里叶系数|ui TyTest|的向量。如在附录1中所述,将X矩阵公式化。MATLAB的内部矩阵乘法和反转应用,用于以求解方程式53。最后,构建由以k为指数的平方差总和SSE组成的向量:
SSEk=(Y-Xβ1sq)T(Y-Xβ1sq) 54
按照Koh′s算法,选择最佳截短指数作为指数,其属于具有最小值的SSE阵列的元素。
在Hansen,P.C.,T.K.Jensen,和G.Rodriguez,An adaptive pruningalgorithm for the discrete L-curve criterion.Journal of Comput and Appl Math,2007.198:第9页中,完整地描述了用于确定截短指数的适当修正算法。总之,所述技术构成离散的L-曲线,它仅仅是解范数
相对于残差范数
的图。Hensen适当修正算法通过在2个阶段中系统地从离散L-曲线除去点,寻找L-曲线的角。在Hansen的正则化工具箱中的角函数执行修正算法,并用于计算截短指数k。Hansen,P.C.,Regularization tools:a MA TLAB package for analysis and solution of discrete ill-posed problems.Numer.Algorithms,1994.6:第35页。
使用所述混合模型的模拟来确定优选的方法,以计算截短指数。对于所有模拟,使用固定流率和体积的造影剂(5ml/s,350mgI/ml造影剂,4秒)作为注射液,模拟在5.1.4.2部分中相同的20个受试者。完成2次模拟运行。第一次使用L-曲线适当修正算法来确定tSVD的最佳截短指数,第二次使用Koh方法来选择截短指数。进行两种方法之间的RMSE对比,以确定优选的方法。
在表15中,给出了20位受试者的3种性能度量RMSE、PDME和EDI的平均值。Koh方法产生预测的增强,具有减小11%的RMS误差、减小22%的最大差异误差和减小10%的平均EDI百分比。并且,当Koh技术计算截短指数时,得到的脉冲响应和重构的测试快速推注曲线更光滑。并且,就该问题而言,Koh方法计算出比适当修正算法更低的值截短指数。基于这些结果,选择Koh方法,用于后续实验。
表11.对比确定tSVD截短指数的2种方法的模拟分析结果
评价方法
一组实验使用混合PBPK模型来产生增强数据,并用于测定在有可加噪声存在下且含有截短的测试快速推注测量向量的非参数估测器性能。接着,使用回顾临床数据(与上述的参数的MLE技术使用的相同的集合)来确定含有人受试者数据的非参数技术的性能。然后,将性能度量与得自上述的参数方法的那些进行对比。
混合模型模拟实验。再次,使用上述的相同方法,使用混合PBPK模型和得自临床数据集的患者人口统计数据,数字地模拟一组20位受试者。仅在左心隔室中的测试快速推注响应用于产生该系统的非参数估测值。20位受试者块具有加入测试快速推注数据中的AWGN(0、0.1、0.25、0.5、1、2、5、10、20HU)。通过用变量长度测试快速推注增强数据进行模拟,研究测试快速推注向量对性能的影响。
临床数据实验。再次使用上述的方法,进行使用回顾临床数据集的实验。使用相同的性能度量RMSE、PDME和EDI,描述估测和预测技术的性能。还进行与使用MLE方法的结果的对比。
使用得自混合模型模拟的模拟数据,使用tSVD来预测对比增强的所有20位受试者的平均RMSE、PDME和EDI结果,显示在图29中。当将相同的诊断注射方案递送给所有20位模拟患者时,在左心隔室的模拟的和预测的增强之间进行对比。这些结果与图21中的那些结果相当,在后者中,MLE方法估测这20位受试者的对比增强。增强值之间的RMSE(平均值±标准差)是7.3±5.1HU,PDME是1.8±1.3%,和EDI是1.6±1.3%。使用含有这些相同数据的MLE,产生了7.78±4.40HU、1.29±1.1%和1.57±1.2%的RMSE、PDME和EDI(参见图21)。
变量长度测试快速推注增强数据(20-40秒)对tSVD估测器性能的影响显示在图30中。所述数据是在每个测试快速推注增强向量长度处取自所有20位受试者的平均值。没有将噪声加入这些实验的测试快速推注增强数据中。
当用AWGN讹误测试快速推注增强数据时,tSVD估测器的性能显示在图31中。使用得自临床数据集的人口统计数据,进行20个模拟。对于所有模拟,在35秒截短测试快速推注向量。将具有0、0.1、0.25、0.5、1.0、2.5、5.0、10和20HU的标准差的可加高斯噪声加入每个20位受试者块的测试快速推注数据中。在该图中一起绘制了tSVD和MLE估测器的平均RMSE值。应当指出,当可加噪声大于2.5HU时,MLE表现得更好。对于低水平的可加噪声(0.1、0.2、0.5和1.0HUσ),tSVD表现得稍微更好。
在图31中给出了tSVD和MLE方法的对比结果,其中使用临床数据作为测试快速推注和诊断性增强数据的源。MLE估测技术分别在14/20、13/20和13/20位受试者中产生了具有更低RMSE、PDME和EDI的增强估测。3种性能度量RMSE、PDME和EDI的总结结果(平均值±SD)是:49.1±24.7HU、14.9±10.5%和13.6±9.6%。
在图33中的箱线图表明,MLE方法在20位受试者中产生更低的RMSE、PDME和EDI,但是差异不是统计上显著的。从箱线图显而易见,使用tSVD的估测增强导致在中位值周围更宽的变异性。在3组度量上的曼-惠特尼U检验揭示了仅PDME结果的显著差异(p<0.05)。
但是,在某些情况下,在2种方法之间证实了所有3种度量的统计显著性。tSVD与MLE相比(p<0.05,Sign-秩检验),受试者5、12和19具有显著更大的平均RMSE和PDME。这些受试者表现出2种数据驱动的估测技术之间的最大偏差。
结果指示,在预测患者特异性的对比增强中,特别当考虑实际的临床实践的约束和实用性时,MLE方法优于tSVD方法。具体地,MLE对测试快速推注增强数据向量长度的变化更稳健的,它具有有利的噪声排斥特征,且它会提供可用于定量患者状态(例如:心输出量估测)的生理学变量的参数估测。如果可得到延伸超过对比增强峰数秒的测试快速推注增强数据,那么tSVD和MLE方法的表现类似。但是,在临床实践中,难以确保在测试快速推注扫描过程中的数据获取常规地含有达到峰值以后的数据点。
然而,tSVD似乎对低水平的AWGN具有优良的噪声免疫,当用大AWGN值讹误测试快速推注数据时,MLE具有更好的噪声免疫。tSVD的主要优点在于,估测仅需要一个TEC曲线,且它在计算方面不太繁重。但是,使用奔腾III(不是特殊硬件)和广泛可得的数字方法,在很短的时间内进行在MLE内的优化。
使用MLE,参数的偏倚和方差是有利的。解空间的可视化也揭示了确定定义的最小值,尽管它们在有些情况下是长的和宽的。因此,有证据表明,当混合模型模拟数据测试MLE方法时,在MLE中使用的代价函数足以基于性能度量来计算患者特异性的对比增强预测所需的参数估测。
在有些情况下,预测和临床数据之间的差一致性可以归于噪音和/或得自临床数据集的不完全测试快速推注TEC。在MLE和tSVD实验中,受试者15的大误差可以直接地归于不完全的测试快速推注曲线。受试者15的TEC绘制在图24中。没有捕获肺动脉和升主动脉的峰增强值。尽管存在这些错失的特征,MLE仍然能够产生对比增强,其形态与临床数据没有过分不相似,尽管在时间上漂移。但是,使用tSVD在该受试者上的3种性能度量都远远大于样品平均值,这说明了有限的测试快速推注数据对tSVD的性能的影响。
MLE结果的第二最差EDI评分是关于受试者7。诊断性增强数据和预测的增强测线的检查表明,增强预测会跟踪诊断性的临床数据,但是偏移40-80HU。该受试者的基线衰减可能低于50HU。
受试者5具有MLE数据的最高RMSE误差。在该受试者的预测的对比增强和临床数据之间存在一般一致性,但是预测的增强测线具有上斜坡,它的开始提前数秒。测试快速推注数据可能是不适当地定时的,或临床数据集扫描延迟被不适当地记录。使用受试者5、12和16观察到含有tSVD预测的最高的RMSE结果。受试者15具有所有3种度量的最差性能,这可归因于该受试者的差测试快速推注数据。
一般而言,在该章中使用的临床扫描数据对于测试识别方法而言不是理想的,因为它们用具有来自测试快速推注方案的不同流率的诊断注射方案产生。如以上所讨论的,注射流率会改变周边隔室流率。在测试快速推注和扫描过程中,通过在一个流率拟合测试快速推注所做出的模型的动力学可以稍微不同于含有不同流率的诊断注射。用于对比识别方法的理想诊断注射方案与混合模型对比所用的类似:所有患者固定的流率和相同的体积。
用于产生诊断性增强数据的方法还必须考虑分析限制,因为TEC不是在构建定时快速推注TEC的相同水平从单水平扫描产生。相反,从造影剂在主动脉中的时空分布,构建出数据。从螺旋数据集产生增强数据,会导致不曾预料的误差和变异性。检查该图中在诊断扫描TEC上的误差棒,表明在TEC数据中存在实质变异性(主要是这些分析的薄片CT数据的选择的结果)。不幸的是,在诊断注射过程中在升主动脉处进行单水平扫描的理想验证情形对于人类而言是不可行的,因为受试者会暴露于过度的辐射和对比。由于这些原因,需要动物和仿真模型来充分地研究造影剂预测模型和对比方案产生技术。动物试验是昂贵的,所以模仿造影剂的运输和分布的经过验证的现实心血管仿真模型具有很大价值。这样的仿真模型描述在,例如,美国专利申请系列号12/397,713。
因为数据驱动的识别估测技术使用从所述系统(使用测试快速推注)导出的数据,它们将产生更大的预测准确度。如预期的,在受试者6、8和12(他们使用如上所述的混合模型具有最大预测误差)中的预测误差下降,因为所述估测器依赖于从所述受试者获取的图像数据。利用更大的临床数据集,可以实现预测误差的统计上显著的减少。
尽管发现MLE方法是优良的方案(例如,当2个TEC可得到时),在临床算法中必须考虑到不可以放置2个ROI的情况,或来自第一个ROI的数据被讹误的情况。因此,预料到,临床算法应当使用tSVD或其它非参数方法作为第二方案,用于在刚刚提及的情况下估测对比增强。例如,可以使用分类来确定TEC数据是否足够MLE技术使用。
总之,2种数据驱动的对比增强方法的对比,揭示了MLE对在临床中遇到的真实考虑更稳健。具体地,当可得到更短的测试快速推注增强数据区段时,MLE可以更好地预测对比增强。但是,tSVD方法仅需要放置一个ROI,且潜在地具有更少的计算负担。最大似然估计被成功地用于鉴别为实现患者特异性的对比增强预测而开发的降阶PBPK模型。使用混合模型和回顾临床数据的比较试验证实,MLE优于tSVD。如以上所讨论的,临床系统、方法或算法可以例如考虑两种数据驱动的技术(参数的和非参数的)用于估测患者特异性的对比增强,因为存在这样的情形:其中不能从测试快速推注产生足够的或适当的数据(例如,2个TEC)来确定参数技术中的参数。在这样的情况下,所述算法可以例如尝试使用tSVD方法来产生对比增强预测。
患者特异性的对比方案/参数产生
上面阐述了用于鉴别和预测人心血管系统中的造影剂药代动力学的患者特异性的和数据驱动的技术。这样的技术不仅适用于预测对比增强,而且适用于提供计算注射方案的方法,所述注射方案会实现为单个患者和操作预期地选择的增强靶标,同时使用最小的造影剂体积。在许多代表性的实施例中,描述了用于计算个体化的注射方案(用于造影剂增强的心胸CT血管造影术)的系统和方法。使用上述的混合PBPK模型作为人数据的替代物,测试了这样的系统和方法。
例如,可以提出问题为非线性的最小化问题,并对输入函数放置非负性约束。这里开发的方案和以前公开的用于计算个体化的对比方案的工作之间的特殊差异是,没有做出努力来迫使预测的对比增强为特定的梯型函数。在许多这样的方案中,导出造影剂注射方案,所述方案尝试在扫描持续时间中实现恒定的、均匀的对比增强。那些方案的一个不利的作用是,与尝试聚焦于实现峰值增强相比,消耗过多的控制能(造影剂体积)来实现均匀的峰值增强。
在过去,CTA扫描的典型获取时间是在20-40秒范围内。由于CT扫描仪的电流产生,CTA数据的获取窗极少超过10秒。由于非常短的扫描获取,设计可最小化或预防倾斜的峰增强测线的注射方案,是并不那么重要的。相反,现代注射方案产生技术可以例如尝试确保,在峰值对比增强过程中发生扫描,并且在数据获取开始时和在扫描获取结束时实现足够的对比增强。在上述对比中使用的临床数据集的扫描持续时间平均为10秒,并在2006年和2007年收集。在过去的几年中,CT技术已经发展,使得那些扫描时间现在比典型扫描持续时间更长。
在图35中,图形地描绘了CTA对比增强的几个目标,和引导本文所述的方案产生方法。混合PBPK模型模拟产生该图所示的对比增强测线。垂直线表示4秒CTA获取窗——当代CT扫描仪的典型扫描持续时间。将扫描定时,使得它位于增强测线上,以确保在扫描获取过程中得到峰值对比浑浊化。在扫描获取(扫描延迟)开始时和在扫描结束时,造影剂衰减大于350HU。如果在该实施例中的扫描持续时间是8秒且希望的增强靶标是300HU,那么注入测线也已经满足目标,假定适当地调节扫描延迟以确保在300HU开始扫描。单个方案产生算法还可以计算单个扫描延迟,因为最佳扫描获取时间作为注射方案、患者的血液动力学和选择的对比增强靶标的函数而变化。
对比增强水平是例如辐射学家可以自由选择的参数,且适当的增强靶标取决于由可疑病理学决定的辐射学家的偏好和任何操作约束。已经出现了一般合意,为了在心胸CTA过程中将脉管系统腔与斑块、血栓和其它病理生理学区分开,250HU的最小对比增强是可接受的。在冠状动脉CTA中,当冠状动脉对比增强大于320HU时,最近证实了用于检测准确的冠状动脉狭窄的增加的灵敏度。在某些情况下,当患者具有肾功能不全或疾病时,可以耐受更低的增强水平,如果这意味着使造影剂体积最小化来帮助预防或减轻肾损伤。
不考虑精确的临床动机,合理的和患者特异性的对比方案产生算法可以例如提供预期地靶向单个患者的对比增强水平的能力。所述算法然后可以尝试在整个扫描持续时间中实现希望的最大对比增强和最小增强。
上面讨论了合理的对比方案产生算法的重要考虑方面。为了实现这样的目标,提出了代价函数,它是对比方案或参数产生方法、系统或算法的许多实施例的组件。通过数字优化,使代价函数最小化。在许多实施例中,最小化的结果是这样的注射方案:其具有足以使成本最小化的最小流率和注射持续时间。由对比方案算法使用的代价函数是
的函数,所述
是通过上述的数据驱动的估测方法产生的左心隔室中的预测的对比增强。它是:
其中TsDly是通过方案产生算法计算的CT研究的扫描延迟,TsDur是CTA获取的扫描持续时间,TPk是最大对比增强的时间,MHU是希望的靶标增强水平。在图35中鉴别了所有参数。将50HU加入峰靶标增强中,认识到对比增强将达到峰值,不做出弄平增强测线的努力。选择50HU作为靶标,因为相差50HU的对比增强的2段主动脉脉管系统通常不是临床上有关的。因为模拟和计算在离散的时间中进行,在代价函数中的时间参数(TsDur、TsDly和TPk)是离散的时间值。
预测的左心增强是参数向量和造影剂注射的函数:
uinj(n)=QinjCinj[u(n)-u(n-Tinj)] (57)
在方程式55中定义的代价函数的最后一个术语会处罚在靶标增强水平以上超过50HU的增强偏差,并被加入代价函数中来预防这样的情况:其中实现了靶HU,但是峰值远远高于必要值(且因此使用比实现靶目标所必需的更多的造影剂)。
放入注射方案产生的背景内,最小化操作被描述为:
其中在向量φ中的辐角是注射方案的流率
和注射的持续时间
它会参数化输入函数方程式57。
最小化是有边界的,但是上下限被放置在注射流率和注射持续时间上。显而易见的下约束是非负性,但是上限是问题特异性的。例如,最大流率可以确定为静脉内导管内径或护理人员偏好的函数。最大注射持续时间是在注射器中可利用的造影剂的最大体积的函数。因为流率和体积的上限由临床约束决定,且是用户可构建的,注射持续时间的上限是最大体积除以最小流率。流率下限设定为3ml/s,最小注射持续时间是8秒。小于3ml/s的注射流率对于CTA而言是不适当的,因为对于大多数患者而言,它们可能导致小于250HU的增强。在3ml/s的最低流率的小于8秒的注射持续时间也可能产生与测试快速推注等同的造影剂体积。
例如,可以在进行数字实验之前研究局部最小值和不连续性在代价函数中的存在,且为此目的,使用以混合PBPK模型产生的模拟对比增强数据,产生由方程式55给出的代价函数的几个表面图。输入函数的参数范围是2-7ml/s(就Qinj而言),以0.1ml/s增量,和6-22秒,以0.5秒增量。使用得自临床数据集的受试者6的人口统计数据来产生对比增强预测yLH(n),在每对注射流率和持续时间。用于产生代价函数空间的图的对比增强靶标和扫描持续时间分别是250HU、350HU靶标和2和8秒。这些值表示典型的上下限,基于电流CT扫描技术和临床偏好。4种情况的三维表面图显示在图36中,解空间的二维投影图呈现为图37中的等值线图。
代价函数表面的检查表明,即使对于1位患者,随着操作参数的变化,解空间的形态学也存在宽泛的变化。8秒扫描显示了2个最小值,但是一个明确定义的全局极小值。对于2秒扫描持续时间情况,全局极小值是独特的,但是它位于长波谷中,这指示数字解算机达到真全局极小值的潜在困难。长波谷指示,可能存在多对可以满足靶标准的注射持续时间和流率(对于等于扫描持续时间的间隔而言,实现希望的靶对比增强)。这些发现会影响下面描述的方案或参数产生方法的开发。
如以上所讨论的,稳健的对比方案会产生适应每位受试者和扫描持续时间的算法,允许辐射学家选择希望的对比增强靶标,并确保这些目标在扫描获取开始和结束时被满足。还应当计算患者特异性的扫描延迟,满足在上下限的参数约束内的要求,并使造影剂的总体积最小化。为了实现这些目标,使用在方程式55和58中定义的代价函数,开发了数字最小化方案,其中含有这样的迭代操作:适应初始条件和参数上边界(当可能时),然后针对靶标增强目标测试解。方案产生技术暗含地使用定时快速推注来产生患者特异性的对比增强估测。所述算法的一般步骤是:(1)给患者注射测试快速推注;(2)使用在
中的初始条件,处理定时快速推注扫描数据,以产生患者特异性的对比增强模型(使用上述技术);和(3)在给定的希望的靶标增强水平(HU)和扫描持续时间下,使用满足方程式58的流率和持续时间计算注射方案。
通过辐射学家的偏好(靶对比增强,MHU)和操作(扫描持续时间,TsDur,最大体积,Vmax和流率,Qmax),确定方程式55的代价函数中的一些参数,同时通过预测的对比增强和输入函数uinj(n),确定最大对比增强的时间(TPk)。
可以考虑T
sDly的参数化,因为它不能由辐射学家或扫描仪操作员容易地定义。它是当扫描获取将要开始时在对比增强
的上斜坡上的时间,且是达到在左心隔室(主动脉)中测得的测试快速推注增强的最大增强
的时间的函数,。
因为测试快速推注具有比完全诊断注射更短的注射持续时间,达到最大测试快速推注对比增强的时间不应当是在诊断扫描过程中开始扫描获取的时间。常规实践是,计算诊断扫描的扫描延迟,作为得自测试快速推注TEC+偏移期TOffset的最大增强时间。在以前的研究中,TOffset是4或6秒,这基于扫描持续时间。对于短扫描,使用更短的偏移时间来尝试将短扫描窗放置在对比增强峰上。当扫描较长时,扫描应当更早地开始,以帮助确保峰和靶阈值拟合在扫描窗内。该扫描定时方案表达为:
其中
是测试快速推注TEC的峰时间。当扫描持续时间超过或等于4秒时,T
Offset是4秒;当扫描持续时间短于4秒时,T
Offset是6秒。算法描述呈现在表16中。
表16方案产生算法
开始A:估测患者模型
1.如果使用MLE技术,使用在第5章中的技术,其基于测试快速推注数据
和
结束A
开始B:计算对比方案
1.获取和储存:MHU、TsDur、VMax、QMax
2.在参数向量上设置上下限
1.φmin=[3 8]
2.
3.初始化TOffset和φ
1.如果TsDur≥4秒
2.TOffset=4秒
3.φ0=[4.5 15]
4.否则
5.TOffset=6秒
6.φ0=[5 13]
7.结束如果
8.m=0
9.
4.当targetMet!=真
1.运行最小化-方程式58直到收敛
i.更新
2.使用估测的患者模型{θHe1D,θHe2D}或hest和uinj(t)(其中 ),计算估测的增强
4.计算TD,从ti时间≥MHU
5.如果TD≥TsDur
i.更新:TsDly=t1,ti的第一个要素
ii.targetMet=真
6.否则
i.增量m
ii.如果m≤5
iv.否则,如果m>5且m≤7
1.如果TOffset==6
a.TOffset=4
2.否则,如果TOffset==4
a.TOffset=2
3.结束如果
v.否则
1.不能满足约束
2.targetMet=真
vi.结束如果
7.结束如果
8.结束While
结束B
VMax是最小化的最大上界约束,并由辐射学家选择。Qmax是通过静脉接入或静脉内导管测量仪确定的单个患者的最大流率。Qmax和Vmax都对最小化常规产生上约束。
通过执行Fleischmann和Hittmair算法,并使用在Fleischmann,D.和K.Hittmair,Mathematical analysis of arterial enhancement andoptimization of bolus geometry for CT angiography using the discrete fouriertransform.J Comput Assist Tomogr,1999.23(3):第474-84页中公开的测试快速推注增强数据作为输入和输出数据,进行与以前公开的结果的对比。表16中的方案产生算法也使用相同的测试快速推注增强数据,并对比用两种方法产生的对比方案。
使用与上述的相同的、用混合模型建立的20个模拟数据集,验证这里提供的方案产生算法。扫描持续时间和靶标增强的多种组合会定义操作数据。从数字最小化和方案算法产生的注射方案用作每位受试者的混合PK模型模拟的输入。如果通过所述算法产生的对比增强测线超过扫描持续时间TsDur的靶标增强(MHU),那么认为结果是成功的。
使用优化工具箱函数fmincon作为数字最优化方法,在MATLAB(R2008b)中实现和执行所述算法。所述fmincon函数执行非线性的有约束的优化,作为连续的二次规划(SQP)问题。SQP技术在解算机的每次迭代处将更大的问题分解为更小的二次规划(QP)子问题。所述方法使用一维搜索技术(在每次迭代时逼近黑塞矩阵,且准牛顿方法会更新拉氏函数)和效益函数来求解每个子问题。解算机的更多细节可以参见:Coleman和Branch,Optimization Toolbox for Use with MA TLAB,User′s Guide,T.Mathworks编,2007。
解算机特异性的约束是:参数的收敛耐受性(TolX)=1E-4,函数的收敛耐受性(TolFun)=1E-4,有限差分子例程的最小增量(DiffMinChange)=0.01,有限差分子例程的最大增量(DiffMaxChange)=2.0,和最大在每次迭代时的函数评价的最大数目是400。迭代的最大数目(maxIter)设定在100。
使用公开的测试快速推注增强数据(如在Fleischmann,D.和K.Hittmair,Mathematical analysis of arterial enhancement and optimization ofbolus geometry ofr CT angiography using the discrete fourier transform.JComput Assist Tomogr,1999.23(3):第474-84页的附录中的数字图表)和希望的增强水平,在MATLAB中执行Fleischmann和Hittmair方法,该方法总结在本文的附录2中,在本说明书的结尾处。得自人受试者的测试快速推注-增强响应(在Fleischmann,D.和K.Hittmair,Mathematical analysis ofarterial enhancement and optimization of bolus geometry for CT angiographyusing the discrete fourier transform.J Comput Assist Tomogr,1999.23(3):第页474-84中报告的增强值)呈现在图38中,并用于产生对比方案。测试快速推注增强响应是由CT扫描仪在降主动脉处测得的短快速推注的典型。将希望的增强定义为200HU水平,其从测试快速推注增强的峰值开始。将0-200HU范围内的递增斜坡和从200HU至0的递减斜坡加入希望的增强测线的前缘和后缘,以产生梯型增强测线。前缘和后缘的斜率与Fleischmann和Hittmair报道的那些匹配。
图39A-D呈现了使用Fleischmann和Hittmair技术计算的注射方案。图39A显示了通过该算法计算的原始方案。所述注射方案在截面上是负的,并超过常用的流率,特别当施用粘稠的(13-20厘泊)、高碘浓度造影剂时。图39B是相同的计算方案,但是具有截短至零的负流率。图39C和D呈现了最大流率设定在7ml/s(任意)的计算方案。图39C和D的方案会截短在下面解释的不同样品点处的输入方案。
达到图38中的测试快速推注增强的峰值增强的时间是14秒,造影剂到达血管(升主动脉)的运输延迟(或造影剂传输时间)是大约8-10秒。就扫描结束之前小于10秒的样品而言,造影剂传输时间会限制注射方案中的造影剂贡献。图39C中的方案在34秒处截短,而图39D的方案在32秒处截短。
使用梯型数值积分来计算每个注射的造影剂总体积,且所述体积注解在图39A-D中。计算这些造影剂体积,以充当与在本论文中呈现的技术进行对比的基础。随着方案的修改,总造影剂体积减少,但是在用于产生这些结果的方案中存在明显的任意性。Fleischmann和Hittmair没有提出优化算法来最小化造影剂体积。
使用从傅里叶解卷积(使用Fleischmann和Hittmair的技术)和图39A-D的计算方案估测的脉冲响应的预测增强测线显示在图40A-D中。就方案A而言,预测的和希望的增强之间的残差最小,就方案D而言最大。
还使用相同的测试快速推注和操作数据,利用在表16中描述的算法,建立了对比方案。因为在Fleischmann,D.和K.Hittmair,Mathematicalanalysis of arterial enhancement and optimization of bolus geometry for CTangiography using the discrete fourier transform.J Comput Assist Tomogr,1999.23(3):第474-84页中仅提供了一条测试快速推注曲线,使用tSVD方法来鉴别患者模型。得到的对比增强曲线(叠合在希望的30秒扫描持续时间上面)显示在图41中。为了实现希望的增强靶标,通过该算法计算的造影剂体积是114.9ml,流率为4.1ml/s。该体积比得自Fleischmann和Hittmair算法的最低体积预测值(图39D)小20ml,且不需要随时间调节流率,所述调节会在外周脉管系统中导致不希望的造影剂快速推注增宽。
使用由混合PBPK模型产生的20个模拟的数据集,在现实条件下测试了方案产生算法。受试者人口统计(身高、体重、性别、年龄)数据与从在上述研究中使用的回顾临床数据集中得到的那些相同。使用在表1和2中定义的生理学参数关联,设定用于模拟运行的混合模型的参数。对于每位受试者,执行所述模型(模拟采样周期=0.01秒/样品),使用在5ml/s的20ml测试快速推注(含有40ml盐水冲洗相)作为输入函数,并将对应的测试快速推注TEC
降低采样。为了更好地模仿临床测试快速推注TEC,在5秒(相对于原点)将向量开窗,并在25秒截短。并且,将零均值AWGN(具有1HU的标准差)加入两种TEC中。
在产生测试快速推注增强曲线以后,将它们用作如上所述开发的MLE算法的输入和输出数据。由在表16中定义的方案产生算法,使用参数化的患者特异性的模型参数(θHe2D)。例如,合乎需要的是,方案产生系统、方法或算法会适应不同的操作特异性的参数和约束。因此,使用模拟的数据来产生在操作设置的组合下的注射方案。用250HU、300HU和350HU的希望靶标增强(MHU)(它们跨被视作心胸CTA可接受的、最低限度地接受的血管对比增强范围),模拟和测试所有20位受试者。在每个增强靶标值,在2、4和8秒的扫描持续时间,进行实验。当用现代MDCT扫描仪(64、128、256、320切片和双源)进行心胸CTA扫描操作时,通常遇到这些扫描持续时间范围。对于所有受试者,最大流率Qmax是7ml/s,最大体积Vmax是150ml。在所有实验中的最小流率是3ml/s,将最小注射持续时间设定为8秒。选择这些值,因为普遍接受的是,低于3ml/s的CTA是不可行的,小于8秒的造影剂快速推注持续时间会导致小于在3ml/s的最小流率的所需快速推注体积(例如21ml)。
通过使用计算的注射方案作为混合PK模型的输入信号,评估方案产生算法的性能。在用计算的注射方案模拟所述模型以后,记录大于希望的MHU的增强的持续时间。将增强持续时间与每个实验集合的扫描持续时间进行对比,将成功定义为在指定的扫描持续时间内比靶标大的增强。因为舍入误差(主要通过TEC向1秒/样品的降低采样引入)和其它随机效应,将可接受性扩展至扫描持续时间±0.5秒。实验构型的流程图呈现在图42中。
来自实验数据集的典型结果的一个实施例显示在图43中。在该情况(得自临床数据集的受试者8)下,扫描持续时间是4秒,MHU是300HU。该算法计算4.52ml/s的诊断注射方案15秒(67.8ml造影剂),并产生超过扫描持续时间的靶标增强的混合模型TEC。
在图44和图45中,显示了在20模拟受试者中表现出超过靶标增强MHU的时间的2个总结图。在两幅图中,水平虚线代表特定实验的扫描持续时间。在图44中,扫描持续时间是2秒,靶标增强是350HU。对于所有20位受试者,增强超过扫描持续时间的MHU(从计算的扫描延迟TsDly开始)。
图45显示了当扫描持续时间是8秒且MHU是250HU时,测试的所有20位受试者的结果。在该实验中,2位受试者具有稍微小于靶标的结果增强时间(受试者3-7.8秒,受试者8-7.6秒),但是在规定的耐受性内。其它受试者具有大于或等于扫描持续时间的结果增强曲线,在某些情况下,具有过度的增强样品时间。但是,与2秒扫描持续时间和350HU靶标结果相比,似乎该优化产生更有效的(最小流率和注射持续时间)方案。所有实验的结果列表显示在表17中,包括计算的诊断注射方案(流率、体积)的总结统计。
表17.得自方案产生实验的总结结果。每行代表在前2列中定义的操作靶标处的混合PK模型的20个模拟。
在所有实验情况下,非线性的求解器发现了可行解,并会聚成解。每位受试者的估测和方案优化的平均执行时间是5.6±1.2秒(奔腾III PC,windows XPPro,2GB RAM)。
使用混合PK模型验证过的技术测试方案产生算法为在不同患者类型和操作设置之间是稳健的。在表17中结果揭示了与对比增强原理相一致的趋势。例如,对于每个靶标增强组,注射持续时间和造影剂体积随着扫描持续时间增加而增加。并且,随着扫描持续时间减少,计算的流率增加,这是意料之中的,因为在短获取下,扫描持续时间应当位于对比增强峰值附近。通过在高流率以短持续时间进行注射,可以实现更大的峰值增强。
在图44所示的结果中,在所有情况下,靶标增强超出几秒。这些结果可以指示,对于短的扫描持续时间,可以例如修改代价函数,以产生更短的注射持续时间。并且,对于短的扫描持续时间,可以重新考虑扫描持续时间在增强曲线上的定位。例如,TsDly可以设置为Tpk减去扫描持续时间的一半。但是,当扫描时间是8秒或更长时,该方案不会产生理想的结果,因为TEC形态学倾向于变得倾斜(向右)更长的注射持续时间。
下面的描述和附图阐述了许多代表性的实施例。当然,本领域技术人员在考虑前述教导以后,会显而易见各种修改、添加和替代设计,且不脱离本发明的范围,所述范围由下述权利要求书而不是前述说明书指示。落入权利要求书的含义和等效范围内的所有变化和变体都被包括在它们的范围内。
附录1
Ostergaard等人应用傅里叶解卷积和其它解卷积方案来提取脑组织的供血动脉和引流静脉之间的“残差”函数或脉冲响应,用于在造影剂增强的(CE)MRI灌注中评估缺血性中风。Ostergaard,L.,等人,Highresolution measurement of cerebral blood flow using intravascular tracer boluspassages.Part I:Mathematical approach and statistical analysis.Magn ResonMed,1996.36(5):第715-25页。鉴于MRI灌注应用不同于当前问题,它与造影剂/患者识别问题具有相似性,因为测量是由使用血管内药剂的生理学血流量系统构成。因此,可以预见到类似的噪声过程。Ostergaard发现,截短的奇异值分解(tSVD)是用于求解组织脉冲响应反转问题的最稳健方法,甚至当SNR<10时。Ostergaard,L.,等人,High resolutionmeasurement of cerebral blood flow using intravascular tracer bolus passages.Part II:Experimental comparison and preliminary results.Magn Reson Med,1996.36(5):第726-36页。最近,Koh等人也证实了tSVD在MR灌注评估中的稳健度,其中使用新颖的技术来选择奇异值分解的截短指数。Koh TS,W.X.,Cheong LH,Lim CCT,Assesment of Perfusion by Dynamic Contrast-Enhanced Imaging Using a Deconvolution Approach Based on Regression andSingular Value Decomposition.IEEE Trans Med Imaging,2004.23(12):第1532-1542页。
tSVD(其它解卷积技术也是如此)的目标是,在给出噪声y测量向量和由输入样品u构成的U矩阵的情况下,估测h。
y=h*u+ε=H·U+ε(A1)
并假定ε是独立同分布的随机过程。在求解h之前,通过奇异值分解,将H分解成单个向量(左边的U和右边的V)和奇异值矩阵(∑):
其中∑是方阵,H的非负奇异值在对角线上。如果将得自方程式A2的H插入方程式A1中,通过代数操作,可以发现h的解。h的一个明显解是:
h=(U∑VT)-1y (A3)
但是,小奇异值可以放大存在于测量向量中的噪声。为了限制噪声的作用,通过忽略大于阈值的奇异值∑,可以减少方程式的秩。h的解向量是:
但是,必须考虑给定问题的截短指数k的适当选择的原理。tSVD是一种求解病态反转问题的正则化技术(非参数方案),其中在数据拟合和解的光滑之间做出折衷。一种标准的正则化技术是Tikhonov的正则化,它求解下述最小化问题(假定线性系统Ax=b的标准符号):
其中λ是正则化参数,其将解范数
相对于残差范数
进行加权。Hansen证实了Tikhonov正则化和tSVD之间的联系,以及如何以SVD的方式表达方程式A5的解:
其中σi是奇异值,λ是正则化参数。总和中的第一个术语称作Tikhonov过滤系数。Hansen,P.C.,The truncated SVD as a method forregularization.BIT,1987.27:第534-55页。如果已知得自tSVD的截短指数k,那么正则化参数是λ=σk,其中σk是在指数k处的奇异值。
L-曲线是正则化参数的不同值的解范数相对于残差范数的重对数坐标图(图47是一个实施例)。使用在L-曲线上的最大曲率的(指数)点,作为解和残差之间的最佳权衡的估测。但是,该技术的一个缺点是,有时角点被噪声遮蔽。
Koh等人提出了一种自动化的技术,其通过使线性曲线逐段拟合SVD的傅里叶系数(|ui Ty|)(相对于奇异值指数)的对数图而发现tSVD的截短指数(并从而发现正则化参数λ=σk)。傅里叶系数的对数图表示称作Picard图。
在有噪声存在下,在Picard图中的傅里叶系数的对数倾向于单调性地下降至这样的点:在该点处,直线的斜率下降,并开始变平(参见图46中的实施例)。偏转指数或转折点则是tSVD的截短指数。在该指数处的奇异值被用于Tikhonov过滤系数中。Koh等人如下确定转折点:使2种不同的线性模型拟合Picard图中的数据,并使用平方差总和标准来确定转折点。
附录2
Fleischmann和Hittmair使用傅里叶方法来执行估测患者/药物传输函数所必需的解卷积。如果ut(n)是离散时间测试快速推注(单位ml/s),yt(n)是代表由扫描仪测得的患者和药物的响应的离散时间信号(亨斯菲尔德单位)。假定造影剂和血管系统表现为线性时间不变系统,通过解卷积发现估测的系统脉冲响应:
其中Z-1代表逆z变换算子。进行离散时间操作,因为每隔1或2秒由CT扫描仪对患者响应信号yt(n)采样。通过将环形目标区(ROI)应用在腹主动脉上面靠近髂动脉处,测量对比增强。
为了减少高频噪声(假定是非生理学的)的影响,Fleischmann和Hittmair通过乘以指数核来过滤输出波谱(m是比例因子):
因为在方程式[sic]的分母中的零或小值可以造成数字不稳定,它们如下确保当Ut(k)是零时H(k)是零:将过滤的输出波谱乘以掩蔽正负号函数sgn(x),后者定义为:
患者/药物系统的估测的传输函数是:
其中:
Tf(k)=sgn(|Ut(k)|) (B6)
掩蔽正负号函数会阻止传输函数估测无限制地增长。
它们的目标是,确定在诊断性CT扫描过程中得到希望的增强响应所需的输入。如下计算理想的造影剂输入函数:
方程式B7的反转的离散傅里叶变换确定及时实现希望的增强水平所需的输入函数。
一般而言,原始的计算的输入函数是不可实现的,因为会产生负流率,且流率经常超过8ml/s。由于注射器、管道和导管的物质强度限制,通过临床和系统约束将流率限制为6-7ml/s。在所述系统中产生的压力是造影剂粘度、导管测量仪(或内径)、管道长度和直径以及患者的血管状态的函数。因此,许多辐射学家根据下述因素减少患者的最大流率:血管接入部位的状况,患者的健康状态,或施用药物的临床医师的信任。为了产生可实现的注入测线,Fleischmann和Hittmair将启发应用于造影剂注射方案,该方案会阻止负流率,并对连续变化的函数四舍五入,使得造影剂注射系统会实现注射。但是,在[9]或[36]中没有描述关于截短的明确算法,尽管在[8]中包括了一些MATHEMATICA代码。