CN101587517A - 用于除去实时pcr荧光数据中的步长不连续的系统和方法 - Google Patents

Abstract

本发明涉及用于除去实时PCR荧光数据中的步长不连续的系统和方法。所述用于自动地除去表示生长过程的数据中的步长不连续的方法，所述方法包括：接收表示生长过程的数据集合；通过将非线性回归处理应用于第一非线性函数以确定第一函数的参数，来计算拟合数据集合的曲线的第一近似；通过将第二回归处理应用于第二非线性函数以确定第二函数的参数，来计算拟合数据集合的曲线的第二近似；为第一和第二近似中的每一个确定信息系数；基于所述信息系数选择所述近似中的一个；针对所选择的近似而确定步长不连续参数的可靠区间；以及如果可靠区间不包括零值，则用对应的步长不连续参数被设置为零的所选择的近似来替换所述数据集合的一部分。

Description

用于除去实时PCR荧光数据中的步长不连续的系统和方法

背景技术

本发明总体上涉及聚合酶链反应系统，以及更具体地涉及用于除去聚合酶链反应数据中的步长不连续的系统和方法。

聚合酶链反应(PCR)是用于酶促地合成或扩增所定义的核酸序列的体外方法。通常，该反应使用两个杂交到相反链并在侧面与模板或要扩增的目标DNA序列相接的低聚核苷酸引子。通过热稳定DNA聚合酶来催化引子的延长。包括模板变性、引子退火和由聚合酶引起的退火引子的延伸的一系列重复周期导致特定DNA片段的指数式积聚。通常，荧光探测器或标记被用于促进扩增处理的检测和量化的处理中。

图1示出了典型的实时PCR曲线，其中相对于用于典型PCR处理的周期数绘出了荧光强度值。在该情形下，在PCR的每个处理周期中对PCR产品的形成进行监视。通常在热循环器中对扩增进行测量，所述热循环器包括用于在扩增反应过程中测量荧光信号的部件和装置。这种热循环器的例子是RocheDiagnotics LightCycler(Cat.No.20110468)。扩增产品例如通过荧光标注的杂交探测器或在某些情形下也通过结合到双链DNA上的荧光染色剂得以被检测，所述荧光标注的杂交探测器在结合到目标核酸上时只发射荧光信号。

用于实时PCR中的荧光染色剂对温度是敏感的。热循环概览图可以包括改变退火温度。退火温度的改变表现为荧光读出值的不连续。图2示出了在周期6之后能看到退火温度的改变的说明性情形。

图3示出了与使用图2的概览图而进行的试验相对应的荧光数据读出。在HEX通道中可以清楚地看到，不连续存在于周期5和6之间。图3取自具有10⁴IU/ml浓度的B19细小-病毒样品。对温度具有更小敏感性的FAM染色剂具有比HEX染色剂更低幅度的跳跃不连续。

但是，如果目标的浓度增加，则不连续在荧光数据中得不到反映是可能的。例如，考虑具有2.9×10¹¹IU/ml浓度的B19细小-病毒样品，即使在周期5之后存在退火温度的改变，在图4a中也看不到不连续。

一些当前的PCR系统实现了用于除去由于退火温度的改变而导致的步长不连续的方法。通常，这样的系统需要针对与退火温度的改变相对应的周期数和所预期的最大荧光改变而分析特定输入参数。这样的方法具有许多限制。首先，它们不能处理高滴度样品。这可以在图4b中看出，在图4b中通过这样的系统对2.9×10¹¹IU/ml浓度下的B19细小样品已经进行了不正确的处理。第二，这样的方法可能只处理不连续点处的小的最大(例如5)的绝对荧光单元变化。已经在多个分析中观察到，不连续处的荧光改变的幅度广泛地不同。建立分析特定的输入以指示最大变化是难以最优化的。第三，这样的方法不能处理基线偏移较高的荧光数据。在这种情形下，所提供的校正可能是不适当的。第四，如果发生尖峰存在于一个周期或小于与退火温度的变化相对应的周期数的距离处，则所提供的校正被已知将导致所报告的结果浓度的变化；即错误的结果。

因此，期望提供用于处理S形(SIGMOID)或生长曲线(尤其是PCR曲线)的、克服上述及其它问题的系统和方法。具体地，所述系统和方法应当以可靠和鲁棒的方式执行温度步长校正。

发明内容

本发明提供用于针对PCR处理过程中可能发生的温度偏移而校正PCR数据的系统和方法。

本发明的实施例消除了实时的聚合酶链反应数据的荧光中的跳跃不连续。不连续的存在可以引起不正确的结果计算，因此建议将其除去。取决于目标的初始复制数，荧光中的不连续在高浓度水平下可能是不存在的。使用八参数非线性模型的可靠区域来准确地估算数据。所述方法克服了之前的用于除去所有浓度水平下的不连续的方法中的困难。

在一个实施例中，通过将非线性回归处理应用于模拟数据集合的非线性函数以确定非线性函数的参数，包括步长不连续参数，来确定拟合所接收的PCR数据集合的曲线的第一近似。非线性函数的一个例子是双S形方程。也通过将回归处理应用于第二非线性函数以确定第二函数的参数，包括步长不连续参数，来确定拟合PCR数据集合的曲线的第二近似。之后，基于为第一近似和第二近似中的每一个所确定的信息系数而选择第一近似和第二近似中的一个。如果针对步长不连续参数而计算的可靠区间包括零值(如可靠区间跨越包括0值的范围(如-1至1))，则不进行步长校正。如果可靠区间不包括零值，则进行步长校正。如果进行步长校正，则步长变化之前的数据曲线部分被替换为所选择的近似的合适部分，以产生偏移校正后的数据集合。在某些方面中，如果近似不能满足拟合标准优度，则不进行步长校正。偏移校正后的数据集合被返回，并且可以被显示或另外用于其他处理。

根据一个实施例，提供一种用于自动地除去表示生长过程的数据中的步长不连续的方法。通常，所述方法包括接收表示生长过程的数据集合，所述数据集合包括多个数据点，每个数据点具有一对坐标值，以及通过将非线性回归过程应用于第一非线性函数以确定第一函数的参数，来计算拟合所述数据集合的曲线的第一近似，所述参数包括步长不连续参数。通常，所述方法还包括：通过将第二回归处理应用于第二非线性函数以确定第二函数的参数，来计算拟合所述数据集合的曲线的第二近似，所述第二函数的参数包括步长不连续参数，为第一近似和第二近似中的每一个确定信息系数，以及基于所述信息系数选择所述近似中的一个。通常，所述方法还包括针对所选择的近似确定步长不连续参数的可靠区间，以及如果可靠区间不包括零值(如可靠区间跨越不包括0值的范围)，则用对应的步长不连续参数被设置为零的所选择的近似代替数据集合的一部分。

根据另一实施例，提供一种计算机可读介质，所述计算机可读介质包括用于控制处理器自动地除去表示生长曲线的数据集合中的步长不连续的代码。通常，所述代码包括用以执行以下处理的指令：接收表示生长过程的数据集合，所述数据集合包括多个数据点，每个数据点具有一对坐标值，以及通过将非线性回归过程应用于第一非线性函数以确定第一函数的参数，来计算拟合所述数据集合的曲线的第一近似，所述参数包括步长不连续参数。通常，所述代码还包括用以执行以下处理的指令：通过将第二回归处理应用于第二非线性函数以确定第二函数的参数，来计算拟合所述数据集合的曲线的第二近似，所述第二函数的参数包括步长不连续参数，为第一近似和第二近似中的每一个确定信息系数，以及基于所述信息系数选择所述近似中的一个。通常，所述代码还包括用以执行以下处理的指令：针对所选择的近似确定步长不连续参数的可靠区间，以及如果可靠区间不包括零值，则用对应的步长不连续参数被设置为零的所选择的近似代替数据集合的一部分。

根据另一实施例，提供动态聚合酶链反应(PCR)系统，通常所述系统包括：动态PCR分析模块，其生成表示动态PCR扩增曲线的PCR数据集合，所述数据集合包括多个数据点，每个数据点具有一对坐标值；以及智能模块。通常，智能模块被适配为：通过计算拟合所述数据集合的曲线的第一近似以及通过计算拟合所述数据集合的曲线的第二近似来处理PCR数据集合，以自动地除去数据集合中的步长不连续，计算拟合所述数据集合的曲线的第一近似是通过将非线性回归处理应用于第一非线性函数以确定第一函数的参数而实现的，所述参数包括步长不连续参数，计算拟合所述数据集合的曲线的第二近似是通过将第二回归处理应用于第二非线性函数以确定第二函数的参数而实现的，所述第二函数的参数包括步长不连续参数。通常，所述智能模块还被适配为：为第一近似和第二近似中的每一个确定信息系数，基于所述信息系数选择所述近似中的一个，针对所选择的近似来确定步长非连续参数的可靠区间，以及，如果可靠区间不包括零值，则用对应的步长不连续参数被设置为零的所选择的近似来代替数据集合的一部分。

参照说明书的剩余部分，包括附图和权利要求，将认识本发明的其它特征和优点。在下面参照附图对本发明的其它特征和优点以及本发明的各个实施例的结构和操作进行了详细的描述。在附图中，相同的附图标记表示相同的或功能相似的元件。

附图说明

图1示出了在PCR处理的情况下扩增曲线的例子。

图2示出了在周期6之后可以看到退火温度的变化的情况。

图3示出了FEM和HEX通道中的、取自具有10⁴IU/ml浓度的B19细小病毒样品的PCR扩增曲线。

图4示出了PCR扩增曲线，该PCR扩增曲线示出：如果目标的浓度增加，则不连续没有反映在荧光数据中是可能的。

图5示出了用于针对PCR扩增曲线中的温度偏移进行校正的处理的一个实施例。

图6A-C示出了根据本发明的一个实施例的用于尖峰识别和替换处理的具体处理流程。

图7示出了包括参数a-g的双S形方程的分解。参数a-g定义双S形曲线的形状和位置。

图8示出了参数(d)对曲线的影响以及(e)的位置对拐点的x值的影响。图8中所有曲线具有相同的参数(除了参数d)。

图9示出了针对不同的参数组的三个曲线形状的例子。

图10示出了根据一方面的用于确定双S形方程的参数(e)和(g)的值的处理。

图11A-B示出了针对初始参数组的Levenberg-Marquardt回归处理的处理流程。

图12示出了来自CAP/CTM HIV监测器的分析的IQS(HEX通道)荧光数据；原始荧光数据如点所示，校正后的数据如实线所示。

图13示出了B19细小病毒目标通道(FAM通道)分析的荧光数据(没有任何实际的目标被扩增)；原始荧光数据如点所示，校正后的数据如实线所示。

图14示出了HBV高纯净(FAM通道)分析的目标的荧光数据；原始荧光数据如点所示，校正后的数据如实线所示。该分析不具有跳跃不连续，所以，不需要对不连续的校正。

图15示出了高浓度(2.9×10¹¹IU/ml)B19细小病毒样品(FAM通道)分析的荧光数据；原始荧光数据如点所示，校正后的数据如实线所示。

图16是示出了可被用于实现本发明的方法和系统的软件资源和硬件资源之间的关系的总体框图的例子。

图17是示出了热循环器装置和计算机系统之间的关系的总体框图的例子。

具体实施例

本发明提供了用于针对在扩增处理过程中可能发生的温度偏移而对PCR扩增曲线和其它生长曲线进行校正的系统和方法。

温度偏移的一个例子是分析过程中某个周期处退火温度的可控变化。通常，该温度偏移发生在由基线区域表示的处理部分期间。该温度变化引起随后的在偏移发生的周期数处的荧光信号的偏移。这里将发生温度偏移的周期称为CAC，即退火变化周期。

在一个实施例中，通过将非线性回归处理应用于模拟数据集合的非线性函数而确定非线性函数的参数(包括步长不连续参数)，来确定拟合PCR数据的曲线的第一近似。非线性函数的一个例子是在下面进行了更具体的描述的双S形方程。还可以通过将回归处理应用于第二非线性函数以确定包括步长不连续参数的第二函数的参数，来确定拟合PCR数据集合的曲线的第二近似。然后，基于为第一和第二近似中的每一个而确定的信息系数来选择第一近似和第二近似中的一个。如果针对步长不连续参数所计算的可靠区间包括零值(如可靠区间跨越包括0值的范围)，则不进行步长校正。如果可靠区间不包括零值，则进行步长校正。如果进行步长校正，则步长变化之前的数据曲线部分被所选择的近似的适当部分代替，以产生偏移校正后的数据集合。在某些方面中，如果近似不满足拟合标准优度，则不进行步长校正。偏移校正后的数据集合被返回，并且可以被显示或另外用于其他处理。例如，偏移校正后的数据集合可以用于确定PCR分析的Ct值。

在PCR处理的情况下的扩增曲线10的一个例子在图1中示出。如图所示，曲线10包括迟滞期区域15和指数期区域25。迟滞期区域15通常被称为基线或基线区域。这样的曲线10包括链接迟滞期区域和指数期区域的受关注的过渡区域20。区域20通常被称为肘弯或肘弯区域。通常，肘弯区域限定了基线的末端和潜在处理的生长或扩增率的过渡。识别区域20中的特定过渡点对于分析潜在处理的行为可能是有用的。在典型的PCR曲线中，识别被称为肘弯值或周期门限(Ct)值的过渡点对于理解PCR处理的效率特性是有用的。

可以提供的类似S形或生长曲线的其它处理包括细菌处理、酶处理和耦合(binding)处理。例如，在细菌生长曲线中，受关注的过渡点被称为迟滞期中的时间λ。根据本发明的产生可以被分析的数据曲线的其它具体处理包括链置换扩增(SDA)处理、基于核酸序列的扩增(NASBA)处理和转录中介扩增(TMA)处理。SDA和NASBA处理和数据曲线的例子可以分别在Wang，Sha-Sha等人的“Homogeneous Real-Time Detection of Signal-Nucleotide Polymorphisms by StrandDisplacement Amp lificafion on the BD ProbeTec ET System”，Clin Chem 200349(10)：1599和Weusten，Jos J.A.M.等人的“Principles of Quantitation of Viral LoadsUsing Nucleic Acid Sequence-based Amplification in Combination WithHomogeneous Detection Using Molecular Beacons”，Nucleic Acids Research，200230(6)：26中找到。所以，虽然本文件的剩余部分将就本发明对PCR曲线的应用性来讨论本发明的实施例和方面，但是应当理解的是，本发明可以应用于与其它处理相关的数据曲线。

如图1所示，典型PCR数据生长曲线的数据可以在二维坐标系统中得以表示，例如，通过PCR周期数定义x轴而积聚的多核苷酸生长指示器定义y轴。通常，如图1所示，积聚生长的指示器是荧光强度值，因为荧光标记物的使用或许是最广泛使用的标记方案。但是，应当理解的是，取决于所使用的特定标记和/或检测方案可以使用其它指示器。丰度或积聚信号生长的其它有用指示器的例子包括发光强度、化学发光强度、生物发光强度、磷光强度、电荷转移、电压、电流、功率、能量、温度、粘性、光散射、放射强度、反射率、透射率和吸光率。周期的定义还可以包括时间、处理周期、单元操作周期和再生周期。

根据本发明，参照图5对用于校正动态PCR扩增曲线的温度偏移的处理100的一个实施例进行描述。在步骤110中，接收或另外获得表示曲线的实验数据集合。还识别发生温度偏移的周期。通常，该周期值为已知的先验值，例如由提供该数据的设备或装置所记录的。图1，3和4示出了绘制的PCR数据集合的例子，其中针对PCR曲线，y轴和x轴分别表示荧光强度和周期数。在某些方面，数据集合应该包括沿着轴线连续的和等间隔的数据。

在处理100在驻留在诸如热循环器的PCR数据采集装置中的智能模块(例如处理器运行指令)中被执行的情况下，在采集数据时可以实时地将数据集合提供给智能模块，或者可以将数据集合存储在存储器中或缓冲器中并在完成试验之后将所述数据集合提供给智能模块。相似地，可以将数据集合经由网络连接(如LAN、VPN、企业内部互联网、国际互联网等)或直接至采集系统的连接(如USB或其它直接有线连接或无线连接)而提供给诸如桌式计算机系统或其它计算机系统的分离的系统，，或将其提供到诸如CD、DVD、软盘、便携式USB驱动等的便携式介质上。在某些方面，数据集合包括具有一对坐标值(或2维向量)的数据点。对于PCR数据，通常一对坐标值表示周期数和荧光强度值。在步骤110中已经接收或获得数据集合之后，可以进一步分析数据集合，用以例如在分析过程中校对温度偏移。

在步骤120中，基于例如从周期1至末周期的数据集合计算第一近似。在该步骤过程中，在一个实施例中，使用具有通过Levenberg-Marquardt(LM)回归处理或其它回归处理而确定的参数的双S形函数来发现表示数据集合的曲线的近似。正如下面更加详细地描述的，在步骤120中，还基于所述数据集合计算不同于第一近似的第二近似。

根据一个实施例，使用其中包含有步长函数的非线性模型的可靠区域。基于双S形模型进行建立(在序列号为11/533,291、名称为“PCR ElbowDetermination Using Curvature Analysis ofa Double Sigmoid”的美国专利申请和序列号为11/861,188、名称为“PCR Elbow Determination Using Quadratic Test forCurvature Analysis of a Double Sigmoid”的美国专利申请中提供了该双S形模型)，根据一个实施例，使用在方程1中指定的函数用作第一非线性模型：

$a+b\·x+\frac{c}{(1+\exp (-d\·(x-e)))(1+\exp (-f\·(x-g)))}+h\·(\mathrm{UnitStep}\left(x\right)-\mathrm{UnitStep}(x-\mathrm{cac}))$

方程1：非线性模型

在方程1中，{a，b，c，d，e，f，g，h}为模型系数，cac为发生退火温度改变的周期数。cac值是已知的，并且是由系统提供的。系数‘h’表示跳跃不连续的幅度。使用非线性回归方法(例如，如在下面将要详细地描述的，Levenberg-Marquardt方法)来推导与非线性模型相对应的系数。通常，非线性回归需要用于合适的整体收敛(其是应用特定的)的初始系数集。在下面更加详细地描述回归处理和参数确定的细节。

在一个实施例中，使用方程2中所定义的UnitStep(单位步长)函数：

$\mathrm{UnitStep}\left(x\right)=\left\{\begin{array}{cc}1& x\&GreaterEqual;0\\ 0& x<0\end{array}\right.$

方程2：UnitStep函数

根据一个实施例，在方程3中指定第二模型：

a+b·x+h·(UnitStep(x)-UnitStep(x-cac))

方程3：附加非线性模型

在一方面中，使用第二非线性模型来准确地表示从没有目标的反应井生成的荧光数据。使用非线性回归方法(例如最小二乘方QRD拟合或其它回归方法)来推导与第二非线性模型相对应的系数。

根据严格的数学观点，并且如从上述实施例的方程中可以看出的，可以将在本发明的方法中所执行的第二函数指定为线性函数，以代替非线性函数。但是，因为在本发明方法中将该第二函数当作非线性函数，所以在本发明申请的情况下将该第二函数指定为非线性函数。

根据一个实施例，在步骤130中，在方程1中所指定的模型和方程3中所指定的模型之间进行选择，在某些方面，使用基于逼近法的非参数信息理论来进行在方程1中所指定的模型和方程3中所指定的模型之间的选择。在一个实施例中，针对每个非线性回归模型(方程1和2)计算如在方程4中所定义的修改的Akaike信息系数(Modified Akaike Information Coefficient，aic)。在该实施例中，假定具有最低系数的回归模型为最佳拟合。在一方面中，从由方程3所指定的模型减去偏差值。该偏差使得能够准确地检测不具有生长的样品。该偏差值可以在从0(无偏差)至大约50的范围内变化。在某些方面中，减去大约为10的偏差。

$\mathrm{aic}=n\&CenterDot;\ln \left(\frac{{\left|\right|\stackrel{\&CenterDot;\&CenterDot;}{y}-y\left|\right|}^2}{n}\right)+2\&CenterDot;m+\frac{2\&CenterDot;m\&CenterDot;(m+1)}{n-m-1}$

方程4：修改的Akaike信息系数

在方程4中，

为预测余差，n是荧光获得周期的数量，m是模型的自由度(例如对于方程1为8，对于方程3为3)。

根据一个实施例，在步骤140中，对拟合优度统计值(例如R²)进行计算以验证模型表达式真正匹配荧光数据。如果该统计值处于预定范围之外，则所述模型不收敛且没有荧光数据的不连续校正被执行。例如，在一方面中，如果R²＜0.8，则没有实现收敛。

该R²值是关于原始曲线(荧光数据)和所估算的模型而进行计算的。在一方面中，给出如下所示的用以计算R²值的方程：

$r^2=1-\frac{{\mathrm{\&Sigma;}}_{i=0}^{n-1}{(f\left(x_i\right)-y_i)}^2}{{\mathrm{\&Sigma;}}_{i=0}^{n-1}{(y_i-\stackrel{\&OverBar;}{y})}^2}$

其中y表示原始数据曲线，f(x)表示所模拟的曲线，y表示原始数据曲线平均值，n表示原始数据曲线的长度(以周期为单位)。

根据一个实施例，在步骤150中，在给定的可靠水平(例如，在90至99.9％之间，诸如95％)处，根据方程1或方程3(取决于选择了哪个)计算h系数的可靠区间。如果该可靠区间包括零值，则模型将系数h作为超定参数。在这种情况下，不识别数据中的不连续且不需要荧光的不连续校正。针对h系数的可靠区间计算取决于所选择的模型和该模型的标准误差。在某些方面中，以如下方程来确定可靠区间的宽度：

$\mathrm{ci}=1.96\&CenterDot;\sqrt{H_h^{-1}\&CenterDot;\mathrm{stdError}},$

其中1.96表示95％处的可靠度，H-1表示针对h系数的Hessian矩阵的逆矩阵，stdError表示所选择的模型的标准误差。以下式来计算Hessian矩阵：

H＝J^T·J

其中J是雅克比(Jacobian)矩阵，J^T是Jacobian转置矩阵。Jacobian是每个荧光周期处的每个模型系数的第一偏导数的估算值。以下式来计算标准误差：

$\mathrm{stdError}=\frac{{\mathrm{\&Sigma;}}_{i=0}^{n-1}{(f\left(x_i\right)-y_i)}^2}{n-2}$

其中y表示原始曲线，f(x)表示所模拟的曲线，n表示曲线的长度。将可靠区间的较低和较高限制计算为h±ci，其中h的值是由估算模型给定的。

应当理解的是，步骤140和150(包括165)是可互换的，且两者都是可选的。

根据一个实施例，在步骤170中，产生步长或偏移校正后的数据集合。在一方面中，通过将系数h设为等于零而使用从所选择的模型获得的值来替换原始荧光强度数据值。在一方面中，为退火温度的变化之前所发生的周期执行替换；这有效地消除了跳跃不连续且确保平滑的过渡。在一方面中，从第一周期开始直到并且包括发生退火温度的变化的周期之后的周期，例如周期1至cac+1，来执行替换。

在步骤175中，返回修改的、偏移校正后的数据集合，例如以用于显示或其他处理。例如，可以对校正后的曲线进行处理以确定Ct值，以及可以将结果(偏移校正后的数据和/或Ct值)返回到例如执行分析的系统，或返回到请求分析的分离的系统。可以利用诸如监视器屏或打印机的显示装置来呈现图形显示，所述显示装置与执行图5的分析的系统相耦合，或者可以将数据提供给用于在显示装置上进行呈现的分离的系统。可以使用校正后的数据集合根据各种方法来确定Ct值。例如，在一方面中，可以使用序列号为11/316315和11/349550的美国专利申请的教导来确定Ct值。

在一个实施例中，提供具有三个不同的开始条件的方程1以导出3组系数，以及提供具有单一组开始条件的方程3以导出单一组系数。在该实施例中，使用四组导出的系数处理方程4以确定具有如上所述的最低aic值的模型。应当认识到的是，不同的开始条件可以用于用来模拟处理的每个非线性方程。在某些方面中，开始条件是凭经验而确定的。

图6A-C的步骤502至524示出了用于近似数据集合的曲线以及确定拟合函数的参数的处理流程。在一个实施例中，使用Levenberg-Marquardt(LM)方法来计算基于曲线点的数据集合的鲁棒曲线近似。LM方法是非线性回归处理；其为使非线性函数和数据集合之间的距离最小化的迭代技术。该处理表现得像是最速下降处理和高斯-牛顿(Gauss-Newton)处理的结合：在当前近似没有很好地拟合时，其表现得像是最速下降处理(较慢但更可靠的收敛)，但是随着当前的近似变得更准确，它则将表现得像是Gauss-Newton处理(较快但较不可靠的收敛)。

通常，LM回归方法包括需要各种输入并提供输出的算法。在一方面中，所述输入包括待处理的数据集合、用于拟合数据的函数以及对所述函数的参数或变量的初始推测。所述输出包括使所述函数和所述数据集合之间的距离最小化的用于所述函数的一组参数。

根据一个实施例，拟合函数为下列形式的双S形：

$f\left(x\right)=a+\mathrm{bx}+\frac{c}{(1+{\exp }^{-d(x-e)})(1+{\exp }^{-f(x-g)})}---\left(5\right)$

选择该方程作为拟合函数是基于其灵活性以及其能能够拟合典型的PCR曲线或其它生长曲线可能采取的不同曲线形状的能力。本领域技术人员应当认识到的是，也可以根据需要而使用上面的拟合函数的变种或其它拟合函数。例如，在一个实施例中使用方程1。

双S形方程(5)具有7个参数：a，b，c，d，e，f和g，而方程1具有8个参数：a，b，c，d，e，f，g和h。可以将该方程分解为常数、倾斜和双S形的总和。双S形自身是两个S形的乘积。图7示出了双S形方程(5)的分解。参数d，e，f和g确定两个S形的形状。为了示出它们对最终曲线的影响，考虑单一S形：

$\frac{1}{1+{\exp }^{-d(x-e)}},---\left(6\right)$

其中参数d确定曲线的“锐度”，参数e确定拐点的x值。图8示出了参数d对曲线的影响以及参数e对拐点的x值的位置的影响。下面，表1描述了参数对双S形曲线的影响。

表1：双S形参数描述

参数	对曲线的影响
		a	x＝0处的y值
b	基线和坪斜(plateau slope)

c	曲线的AFI
		d	第一S形的“锐度”(参见图9)
e	第一S形的拐点的位置(参见图10)
		f	第二S形的“锐度”
g	第二S形的拐点的位置

在一方面中，为了防止曲线具有不切实际的形状，应当约束双S形方程的“锐度”参数d和f。因此，在一方面中，d＜-1或d＞1.1时的任何迭代或者f＜-1或f＞1.1时的任何迭代被认为是失败的。在其它方面中，可以使用对参数d和f的不同约束。

因为Levenberg-Marquardt算法是迭代算法，所以通常需要对待拟合的函数的参数的初始推测。初始推测越好，则近似将越好，以及较不可能发生算法朝向局部最小值收敛。由于双S形函数的复杂性以及PCR曲线或其它生长曲线的各种各样的形状，针对每个参数的一个初始推测可能不足以防止算法有时会朝向局部最小值收敛。因此，在一方面中，输入多(例如三或更多)组初始参数并保存最佳结果。在一方面中，横跨所使用的多组参数，大多数参数保持不变；只有参数c，d和f针对多个参数组的中每一组来说可能是不同的。图9示出了针对不同参数组的三个曲线形状的例子。三组参数的选择指示表示PCR数据的曲线的三个可能的不同形状。应当理解的是，可以对多于三组的参数进行处理并保存最佳结果。

如图6A-C所示，在步骤510中对LM方法的初始输入参数进行识别。可以由操作员输入这些参数，或者可以计算这些参数。根据一方面，根据如下所讨论的步骤502、504和506来确定或设置参数。

初始参数(a)的计算：

参数(a)是基线的高度；其值对于所有的初始参数组来说是相同的。在一方面中，在步骤504中，根据数据集合将参数(a)赋值为第3低的y轴值，例如荧光值。这为鲁棒计算作准备。在其它方面中，当然可以将参数(a)赋值为期望的任何其它荧光值，诸如最低的y轴值、第二低的值等。

初始参数(b)的计算：

参数(b)是基线和坪的倾斜。其值对于所有的初始参数组来说是相同的。在一方面中，在步骤502中，将0.01的静态值赋值给(b)，这是因为理想地不应该有任何倾斜。在其它方面中，可以向参数(b)赋值不同的值，例如，在0至大约0.5的范围内变化的值。在一方面中，值(b)表示从CAC+1至基线末端的基线倾斜。

初始参数(c)的计算：

参数(c)表示坪的高度减去基线的高度，其被表示为绝对的荧光增加，或AFI。在一方面中，对于第一组参数，c＝AFI+2，而对于最后两参数，c＝AFI。这在图9中被示出，其中对于最后两组参数，c＝AFI。对于第一组参数，c＝AFI+2。该变化是由于由第一组参数所模拟的曲线的形状，其不具有坪。

参数(d)和(f)的计算：

参数(d)和(f)定义两个S形的锐度。由于无法基于针对这些参数的曲线给出近似，因此在一方面中，在步骤502中使用三个静态代表值。应当理解的是，其它静态或非静态值可以用于(d)和/或(f)。这些对模拟有关所遇到的PCR曲线的最常见的形状。下面，表2示出了针对如图9所示的不同参数组的(d)和(f)的值。

表2：参数d和f的值

参数组编号	d的值	f的值
			1	0.1	0.7
2	1.0	0.4
			3	0.35	0.25

参数(e)和(g)的计算：

在步骤506中，确定参数(e)和(g)。参数(e)和(g)定义了两个S形的拐点。在一方面中，它们两个跨越所有的初始参数组都采用相同的值。参数(e)和(g)可以具有相同的或不同的值。为了找到近似，在一方面中，使用在诸如荧光的强度的平均值(其不是尖峰)之上的第一个点的x值。图10示出了用于根据该方面确定(e)和(g)的值的处理，且在下面对其进行了讨论。用于根据该方面确定参数(e)和(g)以及其它参数的值的处理的更详细的描述可以在序列号为11/316315的美国专利申请中找到。

参考图10，起初，确定曲线的平均值(例如荧光强度)。接下来，识别在该平均值之上的第一个数据点，然后确定：

a.该点是否不是位于曲线的起始部分的附近，例如前5个周期内；

b.该点是否不是位于曲线的末端附近，例如最后5个周期内；以及

c.该点附近(例如围绕该点半径为2个点之内)的导数是否没有显示出任何符号变化。如果显示出符号变化，则该点可能是尖峰，因此应该抛弃该点。

在另一个实施例中，在确定系数(e)和(g)时只考虑cac周期之后的周期的荧光值。所以，为了找到近似，使用在周期cac至最后一个周期的强度(例如荧光)的平均值之上的第一个点(其不是尖峰)的x值。首先，确定(周期cac至最后一个周期的)曲线的平均值，以及确定平均值之上的第一个点。如果没有找到这样的点，在一方面中则考虑原始曲线长度的三分之二。

下面，表3示出了根据一方面的如图9中所使用的初始参数的例子。

表3：初始参数值：

初始参数组编号	1	2	3
				a的值	第三低的荧光值	第三低的荧光值	第三低的荧光值
b的值	0.01	0.01	0.01
				c的值	第三高的荧光值-a+2	第三高的荧光值-a	第三高的荧光值-a
d的值	0.1	1.0	0.35
				e的值	荧光平均值之上的第一个非尖峰点的X	荧光平均值之上的第一个非尖峰点的X	荧光平均值之上的第一个非尖峰点的X
f的值	0.7	0.4	0.25
				g的值	荧光平均值之上的第一个非尖峰点的X	荧光平均值之上的第一个非尖峰点的X	荧光平均值之上的第一个非尖峰点的X

回到图6A-C，一旦在步骤510中设置了所有的参数，就使用输入的数据集合、函数和参数来执行LM处理520。传统上，Levenberg-Marquardt方法用于解决非线性最小二乘问题。传统的LM方法计算定义为曲线近似和数据集合之间的误差的平方和的距离测量。但是，当最小化平方和时，其给予异常值重要的权重(这是因为异常值的距离大于非尖峰数据点的距离)，从而常常导致不适当的曲线或较不理想的曲线。因此根据本发明的一方面，通过最小化绝对误差的总和来计算近似和数据集合之间的距离，这是因为这样不会给予异常值如此多的权重。在该方面中，由下式给出近似和数据之间的距离：

distance＝∑|y_data-y_{approximation}|    (7)

如上，在一方面中，如步骤522和524中所示，对多(例如3)组初始参数中的每组进行输入和处理，并保存最佳结果，其中最佳结果是在方程(7)中提供最小或最少距离的参数组。在一方面中，跨越多组参数，大多数参数保持不变；对于每组参数只有c，d和f可能是不同的。应当理解的是，可以使用任意数量的初始参数组。

图11A-B示出了根据本发明的针对参数组的LM处理520的处理流程。如上面所解释的，Levenberg-Marquardt方法可以表现得像是最速下降处理或像是Gauss-Newton处理。其行为取决于阻尼因子λ。λ越大，Levenberg-Marquardt算法就表现得越像最速下降处理。另一方面，λ越小，Levenberg-Marquardt算法就表现得越像Gauss-Newton处理。在一方面中，λ是从0.001开始的。应当理解的是，λ可以从任何其他值处开始，诸如从约0.000001至约1.0。

如前所述，Levenberg-Marquardt方法是迭代技术。根据一方面，如图11A-B所示，在每次迭代过程中进行下面的处理：

1.计算之前近似的Hessian矩阵(H)。

2.计算之前近似的转置Jacobian矩阵(J^T)。

3.计算之前近似的距离向量(d)。

4.以当前的阻尼因子λ来扩增Hessian矩阵的对角线：

H_aug＝Hλ    (8)

5.对扩增后的方程求解：

H_augx＝J^Td    (9)

6.将扩增后的方程的解x加到函数的参数上。

7.计算新的近似和曲线之间的距离。

8.如果在该新的参数组情况下的距离小于在先前参数组情况下的距离，则：

●认为迭代是成功的。

●保存或存储该新的参数组。

●减小阻尼因子λ，例如以因子10来减小。

如果在该新的参数组情况下的距离大于在先前参数组情况下的距离，则：

●认为迭代是不成功的。

●丢弃该新的参数组。

●增大阻尼因子λ，例如以因子10来增大。

在一方面中，图11A-B的LM处理进行迭代，直到达到下面的标准中的一个为止：

1.其已经运行达到了指定的迭代次数N。该第一标准防止算法无穷地迭代。例如，在一方面中，如图10所示，默认的迭代值N为100。对于算法来说100次迭代应该是足以收敛的(如果它可以收敛的话)。一般地，N可以在从小于10至100或更大的范围内变化。

2.两次成功的迭代之间的距离的差小于阈值，例如0.0001。当该差变得很小时，已经达到了期望的精度，并且继续进行迭代是无意义的，因为解不会变得显著地更好。

3.阻尼因子λ超过指定值，例如大于10²⁰。当λ变得很大时，算法不会比当前的解更好地收敛，因此继续进行迭代是无意义的。一般地，该指定值可以显著地小于或大于10²⁰。

根据另一实施例，所使用的Levenberg-Marquardt方法等同于在下面的方程10中所指定的表达式。假设向量函数具有初始条件z₀和精度水平tol(例如默认设置为10^-3)，方程10有效。函数J(z)被定义为函数f(z)的Jacobian矩阵。变量maxiter为最大的迭代次数(例如默认设置为100)。变量λ在回归中是阻尼因子。

方程10：Levenberg-Marquardt非线性回归

利用QR分解进行方程10的倒置。在一方面中，QR分解类似于Anderson的DGELS方法，E.(ed)，“LAPACK Users’Guide”第三版。Philadelphia，PA：Society for Industrial and Applied Mathematics，1999。如果不能确定倒置，则用2(或其它值)乘以λ，并且回归继续。

在本发明的方法、系统和计算机可读介质的某些实施例中，为第一近似和第二近似中的每一个确定信息系数可以包括为第一近似和第二近似中的每一个计算Akaike信息系数(AIC)值。

在某些实施例中，根据本发明的方法、系统和计算机可读介质还可以包括从针对第二非线性函数所计算的AIC值中减去偏差值。

在本发明的方法、系统和计算机可读介质的某些实施例中，基于信息系数选择近似中的一个可以包括选择具有最低AIC值的近似。

在本发明的方法和计算机可读介质的某些实施例中，数据集合可以表示动态聚合酶链反应(PCR)处理的生长曲线，其中坐标值对可以表示扩增后的多核苷酸的积聚和周期数。

在本发明的方法和计算机可读介质的某些实施例中，扩增后的多核苷酸的积聚可以由荧光强度值、发光强度值、化学发光强度值、磷光强度值、电荷转移值、生物发光强度值或吸光率值中的一个来表示。

在本发明的方法和计算机可读介质的某些实施例中，数据集合可以表示动态聚合酶链反应(PCR)处理、细菌处理、酶处理或耦合处理的生长曲线。

在根据本发明的系统的某些实施例中，所替换的数据集合部分可以包括从第一数据点开始直到且包括步长不连续发生处的数据点的数据集合部分。

另外，智能模块还可以被适配为在选择之后针对所选择的近似计算拟合值的优度，并且仅当拟合值的优度超过阈值时才继续。

在根据本发明的系统的某些实施例中，智能模块还可以被适配为从针对第二非线性函数所计算的AIC值中减去偏差值。

在根据本发明的系统的某些实施例中，基于信息系数选择近似中的一个可以包括选择具有最低AIC值的近似。

在根据本发明的系统的某些实施例中，动态PCR分析模块可以驻留在动态热循环器装置中，其中智能模块可以包括可通信地耦合到分析模块的处理器。

在根据本发明的系统的某些实施例中，智能模块可以包括驻留在通过网络连接和直接连接中的一个而耦合到分析模块的计算机系统中的处理器。

在根据本发明的系统的某些实施例中，所替换的数据集合部分可以包括在步长不连续发生处的数据点之后的数据点。

在根据本发明的计算机可读介质的某些实施例中，所替换的数据集合部分可以包括在步长不连续发生处的数据点之后的数据点。

在根据本发明的方法的某些实施例中，所替换的数据集合部分可以包括在步长不连续发生处的数据点之后的数据点。

在某个实施例中，可以通过使用传统的个人计算机系统来实现根据本发明的方法，所述传统的个人计算机系统包括但并不局限于：用以输入数据集合的输入装置，诸如键盘、鼠标等；用以表现曲线的区域中的特定关注点的显示装置，诸如监视器；执行方法中的每个步骤所必需的处理装置，诸如CPU；诸如调制解调器的网络接口、用以存储数据集合的数据存储装置、在处理器上运行的计算机代码等。此外，还可以在根据本发明的PCR装置中实现所述方法。

根据本发明的系统的例子在图16-17中被显示。图16示出了解释可以用于实现本发明的方法和系统的软件资源和硬件资源之间的关系的整体框图。在图17中示出的系统包括可以位于热循环器装置中的动态PCR分析模块和作为计算机系统的一部分的智能模块。将数据集合(PCR数据集合)从分析模块经由网络连接或直接连接传送至智能模块，或反之亦然。例如，可以根据图5、6A-B、10和11A-B所示出的流程图来处理数据集合。例如根据图16所描述的流程图，可以通过存储在计算机系统的硬件上的软件来方便地实现这些流程图。参照图16，计算机系统(200)可以包括例如用于接收在PCR反应过程中获得的荧光数据的接收装置(210)、用于根据本发明的方法处理所述数据的计算装置(220)、用于根据由计算装置获得的结果来替换所述数据的一部分的应用装置(230)和用于将结果显示在计算机屏幕上的显示装置(240)。图17示出了热循环器装置和计算机系统之间的交互作用。该系统包括可以位于热循环器装置中的动态PCR分析模块和作为计算机系统的一部分的智能模块。将数据集合(PCR数据集合)经由网络连接或直接连接从分析模块传送至智能模块，或反之亦然。可以通过在处理器上运行且存储在智能模块的存储装置上的计算机代码根据图16来处理数据集合，并且在处理之后可以将所述数据集合传送回分析模块的存储装置，其中可以将修改后的数据显示在显示装置上。

如上面所解释的，本发明的系统和方法对于除去聚合酶链反应数据中的步长不连续来说是有用的。例如，当使用荧光数据来监视聚合酶链反应时，本发明的系统和方法提供更准确的数据。这样的数据不仅对监视反应来说是有用的，而且还提供诸如在PCR过程中对扩增的目标核苷酸的量化或根据所获得数据来适配PCR的反应条件的技术效果。

提供下面的例子和图形以帮助对本发明的理解，本发明的真实的范围在所附权利要求中进行了阐述。应当理解的是，在不超出本发明的精神的情况下可以对所阐述的程序作出修改。

例子

例1

这是来自CAP/CTM HIV监视器分析的IQS(HEX通道)荧光数据。该分析图表在周期15之后具有跳跃不连续。图12示出了原始荧光数据(点)和校正后的数据(实线)。下面，表4示出了所估算的系数和对应的可靠区间。例1的R²大于0.99。

表4：例1的系数图表

例2

这是在没有任何实际的目标被扩增的情况下的B19细小-病毒目标通道(FAM通道)。该分析图表在周期5之后具有跳跃不连续。图13示出了原始荧光数据(点)和校正后的数据(实线)。下面，表5示出了所估算的模型系数和它们的对应可靠区间。例2的R²大于0.94。

表5：例2的系数图表

例3

这是HBV高纯净的目标(FAM通道)。该分析不具有跳跃不连续。在该例子中系数h的可靠区间将包括零。所以，不需要不连续的校正。图14示出了原始数据(点)和校正后的数据(实线)。下面，表6示出了所估算的模型系数和它们的对应可靠区间。例3的R²大于0.99。

表6：例3的系数图表

例4

这是高浓度(2.9×10¹¹IU/ml)B19细小-病毒样品(FAM通道)。虽然分析图表在周期5之后具有跳跃不连续，但是由于高复制数，看不到该不连续。上面所述的方法能够很好地处理该例子。图15示出了原始数据(点)和校正后的数据(实线)。在该例子中系数h的可靠区间包括零。下面，表7示出了所估算的模型系数和它们的对应可靠区间。例4的R²大于0.99。

表7：例4的系数图表

虽然已经通过例子并就具体实施例对本发明进行了描述，但是应当理解的是，本发明并不局限于所公开的实施例。

Claims (15)

1.一种用于自动地除去表示生长过程的数据中的步长不连续的方法，所述方法包括：

接收表示生长过程的数据集合，所述数据集合包括多个数据点，每个数据点具有一对坐标值；

通过将非线性回归处理应用于第一非线性函数以确定第一函数的参数，来计算拟合所述数据集合的曲线的第一近似，所述参数包括步长不连续参数；

通过将第二回归处理应用于第二非线性函数以确定第二函数的参数，来计算拟合所述数据集合的曲线的第二近似，所述第二函数的参数包括步长不连续参数；

为所述第一近似和第二近似中的每一个确定信息系数；

基于所述信息系数选择所述近似中的一个；

针对所选择的近似而确定所述步长不连续参数的可靠区间；以及

如果所述可靠区间不包括零值，则用对应的步长不连续参数被设置为零的所选择的近似来替换所述数据集合的一部分。

2.如权利要求1所述的方法，其中第一非线性回归处理是Levenberg-Marquardt(LM)回归处理，其中第一非线性函数是双S形函数。

3.如权利要求2所述的方法，其中所述双S形函数是如下形式：

$a+b\&CenterDot;x+\frac{c}{(1+\exp (-d\&CenterDot;(x-e)))(1+\exp (-f\&CenterDot;(x-g)))}+h\&CenterDot;(\mathrm{UnitStep}\left(x\right)-\mathrm{UnitStep}(x-\mathrm{cac}))$

其中UnitStep是如下形式：

$\mathrm{UnitStep}\left(x\right)=\left\{\begin{array}{cc}1& x\&GreaterEqual;0\\ 0& x<0\end{array}\right.,$ 以及

其中计算第一近似包括迭代地确定函数的参数a、b、c、d、e、f、g和h中的一个或多个。

4.如权利要求1所述的方法，其中第二非线性函数具有如下形式：

a+b·x+h·(UnitStep(x)-UnitStep(x-cac))，其中UnitStep是如下形式：

$\mathrm{UnitStep}\left(x\right)=\left\{\begin{array}{cc}1& x\&GreaterEqual;0\\ 0& x<0\end{array}\right.,$ 以及其中，x是周期数，cac是不连续发生处的周期，其中a、b和h是所述参数，h是所述步长不连续参数。

5.如权利要求1所述的方法，其中所述数据集合的被替换的部分包括从第一个数据点开始直到且包括步长不连续发生处的数据点的数据集合部分。

6.如权利要求1所述的方法，还包括在选择步骤之后：

针对所选择的近似而计算拟合值的优度，以及仅当拟合值的优度超过阈值时才继续。

7.一种计算机可读介质，其包括用于控制处理器自动地除去表示生长曲线的数据集合中的步长不连续的代码，所述代码包括用以执行如下处理的指令：

接收表示生长过程的数据集合，所述数据集合包括多个数据点，每个数据点具有一对坐标值；

通过将非线性回归处理应用于第一非线性函数以确定第一函数的参数，来计算拟合所述数据集合的曲线的第一近似，所述参数包括步长不连续参数；

通过将第二回归处理应用于第二非线性函数以确定第二函数的参数，来计算拟合所述数据集合的曲线的第二近似，所述第二函数的参数包括步长不连续参数；

为所述第一近似和第二近似中的每一个确定信息系数；

基于所述信息系数选择所述近似中的一个；

针对所选择的近似而确定所述步长不连续参数的可靠区间；以及

如果所述可靠区间不包括零值，则用对应的步长不连续参数被设置为零的所选择的近似来替换所述数据集合的一部分。

8.如权利要求7所述的计算机可读介质，其中第一非线性回归处理是Levenberg-Marquardt(LM)回归处理，其中第一非线性函数是双S形函数。

9.如权利要求8所述的计算机可读介质，其中所述双S形函数是如下形式：

$a+b\&CenterDot;x+\frac{c}{(1+\exp (-d\&CenterDot;(x-e)))(1+\exp (-f\&CenterDot;(x-g)))}+h\&CenterDot;(\mathrm{UnitStep}\left(x\right)-\mathrm{UnitStep}(x-\mathrm{cac}))$

其中UnitStep是如下形式：

$\mathrm{UnitStep}\left(x\right)=\left\{\begin{array}{cc}1& x\&GreaterEqual;0\\ 0& x<0\end{array}\right.,$ 以及

其中计算第一近似包括迭代地确定函数的参数a、b、c、d、e、f、g和h中的一个或多个。

10.如权利要求7所述的计算机可读介质，其中所述第二非线性函数具有如下形式：

a+b·x+h·(UnitStep(x)-UnitStep(x-cac))，其中UnitStep是如下形式：

$\mathrm{UnitStep}\left(x\right)=\left\{\begin{array}{cc}1& x\&GreaterEqual;0\\ 0& x<0\end{array}\right.,$ 以及其中x是周期数，cac是不连续发生处的周期，其中a、b和h是所述参数，h是所述步长不连续参数。

11.如权利要求7所述的计算机可读介质，其中所述数据集合的被替换的部分包括从第一个数据点开始直到且包括步长不连续发生处的数据点的数据集合部分。

12.一种动态聚合酶链反应(PCR)系统，包括：

动态PCR分析模块，其生成表示动态PCR扩增曲线的PCR数据集合，所述数据集合包括多个数据点，每个数据点具有一对坐标值；以及

智能模块，其被适配为通过以下操作来处理所述PCR数据集合以自动地除去所述数据集合中的步长不连续：

通过将非线性回归处理应用于第一非线性函数以确定第一函数的参数，来计算拟合所述数据集合的曲线的第一近似，所述参数包括步长不连续参数；

通过将第二回归处理应用于第二非线性函数以确定第二函数的参数，来计算拟合所述数据集合的曲线的第二近似，所述第二函数的参数包括步长不连续参数；

为所述第一近似和第二近似中的每一个确定信息系数；

基于所述信息系数选择所述近似中的一个；

针对所选择的近似而确定所述步长不连续参数的可靠区间；以及

如果所述可靠区间不包括零值，则用对应的步长不连续参数被设置为零的所选择的近似来替换所述数据集合的一部分。

13.如权利要求12所述的系统，其中所述第一非线性回归处理是Levenberg-Marquardt(LM)回归处理，其中第一非线性函数是双S形函数。

14.如权利要求13所述的系统，其中所述双S形函数是如下形式：

$a+b\&CenterDot;x+\frac{c}{(1+\exp (-d\&CenterDot;(x-e)))(1+\exp (-f\&CenterDot;(x-g)))}+h\&CenterDot;(\mathrm{UnitStep}\left(x\right)-\mathrm{UnitStep}(x-\mathrm{cac})),$

其中UnitStep是如下形式：

$\mathrm{UnitStep}\left(x\right)=\left\{\begin{array}{cc}1& x\&GreaterEqual;0\\ 0& x<0\end{array}\right.,$ 以及

其中计算第一近似包括迭代地确定函数的参数a、b、c、d、e、f、g和h中的一个或多个。

15.如权利要求12所述的系统，其中所述第二非线性函数具有如下形式：a+b·x+h·(UnitStep(x)-UnitStep(x-cac))，

其中UnitStep是如下形式：

$\mathrm{UnitStep}\left(x\right)=\left\{\begin{array}{cc}1& x\&GreaterEqual;0\\ 0& x<0\end{array}\right.,$ 以及其中x是周期数，cac是不连续发生处的周期，其中a、b和h是所述参数，h是所述步长不连续参数。

Images