Implantable pacemakers control and optimization via fractional calculus approaches
View Patent ↗Method and system for non-linear modeling of physiological behavior, such as R-R intervals, in implantable devices, such as a rate responsive pacemakers, comprising a comprehensive modeling and optimization methodology based on fractional calculus and constrained finite horizon optimal control theory that allows for allows for fine-grain optimization of pacemaker response to heart rate variations; and the theoretical basis on which a hardware implementation of the fractional optimal controller that can respond to changes in the heart rate dynamics. Present invention describes a fractal approach to pacemaker control based on the constrained finite horizon optimal control problem. This is achieved by modeling the heart rate dynamics via fractional differential equations. Also, by using calculus of variations, the invention describes how the constrained finite horizon optimal control problem can be reduced to solving a linear system of equations. Finally, the invention describes the theoretical basis on which a hardware implementation become possible.
1. A method for non-linear (fractional dynamics) modeling of physiological behavior measured by an implantable device, comprising the steps of:
measuring the physiological behavior;
determining magnitude of deviation between the measured physiological behavior and a reference value y ref ;
modeling dynamics of the measured physiological behavior with continuous time fractional differential equations to identify parameters of a non-linear fractal model;
selecting a constrained fractal optimal control problem;
deriving optimality conditions for the constrained fractal optimal control problem when the magnitude of deviation exceeds a predetermined threshold;
discretizing the optimality conditions for the constrained fractal optimal control problem; and
solving the constrained fractal optimal control problem corresponding to the optimality conditions to identify an optimal physiological behavior to operate the implantable device.
2. The method according to claim 1 , wherein in the physiological behavior is an R-R interval.
3. The method according to claim 1 , wherein the implantable device is a pacemaker.
4. The method according to claim 1 , wherein the step of deriving optimality conditions for the constrained fractal optimal control problem comprises using wavelets and linear regression.
5. The method according to claim 1 , wherein the step of selecting a constrained fractal optimal control problem comprises the step of selecting a finite time fractal optimal control with integral of squared tracking error criterion.
6. The method according to claim 1 , wherein the step of selecting a constrained fractal optimal control problem comprises the step of selecting a finite time fractal optimal control with integral of squared time multiplied by squared tracking error criterion.
7. The method according to claim 1 , wherein the step of determining magnitude of deviation between the measured physiological behavior and a reference value comprises analyzing the magnitude of a state variable y(t) minus the reference valve y ref .
8. The method according to claim 1 , wherein the state variable y(t) is an R-R interval.
9. The method according to claim 1 , wherein the step of deriving optimality conditions for the constrained fractal optimal control problem when the magnitude of deviation exceeds a predetermined threshold, comprises the step of solving the following equations:
min
f
(
t
)
∫
t
i
t
f
{
w
2
[
y
(
t
)
-
y
ref
(
t
)
]
2
+
z
2
f
2
(
t
)
}
ⅆ
t
ⅆ
α
y
(
t
)
ⅆ
t
α
=
a
(
t
)
y
(
t
)
+
b
(
t
)
f
(
t
)
,
y
(
t
i
)
=
y
0
,
y
(
t
f
)
=
y
0
,
0
≤
y
min
≤
y
(
t
)
≤
max
≤
1
f
min
≤
f
(
t
)
f
max
where y(t) represents the heart rate activity seen as a state variable, y ref (t) denotes the reference value, f(t) is the pacing frequency, w and z are the weighting coefficients for the quadratic error and magnitude of the control signal, respectively, in the cost function, α is the exponent of the fractional order derivative characterizing the dynamics of the heart rate activity y(t), a(t) and b(t) are weighting coefficients for the heart activity and pacing frequency, y min and y max are the minimum and maximum bounds on heart rate activity y(t), y(t i ) is the initial condition, y(t f ) is the final condition, f min and f max are the minimum and maximum allowed bounds on pacing frequency f(t);
L
=
∫
t
i
t
f
{
w
[
y
(
t
)
-
y
ref
(
t
)
]
2
2
+
z
f
2
(
t
)
2
+
β
1
(
f
-
f
min
-
ξ
1
)
--
λ
[
ⅆ
α
y
(
t
)
ⅆ
t
α
-
a
(
t
)
y
(
t
)
-
b
(
t
)
f
(
t
)
]
+
β
2
(
f
max
-
f
-
ξ
2
)
}
ⅆ
t
where λ, β 1 , and β 2 are the Lagrange multipliers associated with the dynamical state equation for y(t) and the constraints imposed on the control variable f, ξ 1 and ξ 2 are the slack variables needed to transform the inequality bounds into equality constraints on the control variable f; and
∂
L
∂
y
+
t
D
t
f
α
∂
L
∂
t
i
D
t
α
y
=
0
,
∂
L
∂
f
=
0
,
∂
L
∂
λ
=
0
,
∂
L
∂
β
1
=
0
,
∂
L
∂
β
2
=
0
where t i D t α and t D t f α are the fractional derivatives operating backward and forward in time, respectively.
10. The method according to claim 1 , wherein the step of solving the constrained fractal optimal control problem corresponding to the optimality conditions to identify an optimal physiological behavior to operate the implantable device, comprises the step of solving the following equations:
y ( t i )= y 0 , y ( t f )= y 0 , 0≦ y min ≦y ( t )≦ y max ≦1
f min ≦f ( t )≦ f max
where y(t) represents the heart rate activity seen as a state variable, f(t) is the pacing frequency, y min and y max are the minimum and maximum bounds on heart rate activity y(t), y(t i ) is the initial condition, y(t f ) is the final condition, f min and f max are the minimum and maximum allowed bounds on pacing frequency f(t); and
∑
i
=
0
k
(
-
1
)
i
Δ
t
α
(
α
i
)
y
(
(
k
-
i
)
Δ
t
)
-
a
(
k
Δ
t
)
y
(
k
Δ
t
)
+
b
(
k
Δ
t
)
z
(
k
Δ
t
)
λ
(
k
Δ
t
)
+
b
(
k
Δ
t
)
z
(
k
Δ
t
)
β
1
(
k
Δ
t
)
-
b
(
k
Δ
t
)
z
(
k
Δ
t
)
β
2
(
k
Δ
t
)
=
0
,
k
=
1
,
…
,
N
β
2
(
k
Δ
t
)
-
β
1
(
k
Δ
t
)
-
b
(
k
Δ
t
)
λ
(
k
Δ
t
)
z
(
k
Δ
t
)
-
ξ
1
(
k
Δ
t
)
=
f
min
,
k
=
1
,
…
,
N
β
2
(
k
Δ
t
)
-
β
1
(
k
Δ
t
)
-
b
(
k
Δ
t
)
λ
(
k
Δ
t
)
z
(
k
Δ
t
)
+
ξ
2
(
k
Δ
t
)
=
f
max
,
k
=
1
,
…
,
N
∑
i
=
0
N
-
k
(
α
i
)
(
-
1
)
i
λ
(
(
k
+
i
)
Δ
t
)
Δ
t
α
-
w
(
k
Δ
t
)
[
y
(
k
Δ
t
)
-
y
ref
(
k
Δ
t
)
]
-
a
(
k
Δ
t
)
λ
(
k
Δ
t
)
+
λ
(
Δ
t
N
)
(
t
f
-
t
i
-
k
Δ
t
)
-
α
Γ
(
1
-
α
)
=
0
,
k
=
N
-
1
,
…
,
0
which computes the y and λ values for a predefined set of discretization steps.
11. The method according to claim 1 , wherein the step of deriving optimality conditions for the constrained fractal optimal control problem when the magnitude of deviation exceeds a predetermined threshold, comprises the step of solving the following equations:
min
f
(
t
)
∫
t
i
t
f
{
w
2
[
t
(
y
(
t
)
-
y
ref
(
t
)
)
]
2
+
z
2
f
2
(
t
)
}
ⅆ
t
ⅆ
α
y
(
t
)
ⅆ
t
α
=
a
(
t
)
y
(
t
)
+
b
(
t
)
f
(
t
)
,
y
(
t
i
)
=
y
0
,
y
(
t
f
)
=
y
0
,
0
≤
y
min
≤
y
(
t
)
≤
max
≤
1
f
min
≤
f
(
t
)
f
max
where y(t) represents the heart rate activity seen as a state variable, y ref (t) denotes the reference values that need to be achieved in terms of heart rate activity, f(t) is the pacing frequency, w and z are the weighting coefficients for the quadratic error and magnitude of the control signal, respectively, in the cost function, α is the exponent of the fractional order derivative characterizing the dynamics of the heart rate activity y(t), a(t) and b(t) are weighting coefficients for the heart activity and pacing frequency, y min and y max are the minimum and maximum bounds on heart rate activity y(t), y(t i ) is the initial condition, y(t f ) is the final condition, f min and f max are the minimum and maximum allowed bounds on pacing frequency f(t);
L
=
∫
t
i
t
f
{
w
[
y
(
t
)
-
y
ref
(
t
)
]
2
2
+
z
f
2
(
t
)
2
+
β
1
(
f
-
f
min
-
ξ
1
)
--
λ
[
ⅆ
α
y
(
t
)
ⅆ
t
α
-
a
(
t
)
y
(
t
)
-
b
(
t
)
f
(
t
)
]
+
β
2
(
f
max
-
f
-
ξ
2
)
}
ⅆ
t
where λ, β 1 , and β 2 are the Lagrange multipliers associated with the dynamical state equation for y(t) and the constraints imposed on the control variable f, ξ 1 and ξ 2 are the slack variables needed to transform the inequality bounds into equality constraints on the control variable f; and
∂
L
∂
y
+
t
D
t
f
α
∂
L
∂
t
i
D
t
α
y
=
0
,
∂
L
∂
f
=
0
,
∂
L
∂
λ
=
0
,
∂
L
∂
β
1
=
0
,
∂
L
∂
β
2
=
0
where t i D t α and t D t f α are the fractional derivatives operating backward and forward in time, respectively.
12. The method according to claim 1 , wherein the step of solving the constrained fractal optimal control problem corresponding to the optimality conditions to identify an optimal physiological behavior to operate the implantable device, comprises the step of solving the following equations:
y ( t i )= y 0 , y ( t f )= y 0 , 0≦ y min ≦y ( t )≦ y max ≦1
f min ≦f ( t )≦ f max
where y(t) represents the heart rate activity seen as a state variable, y ref (t) denotes the reference values that need to be achieved in terms of heart rate activity, f(t) is the pacing frequency, w and z are the weighting coefficients for the quadratic error and magnitude of the control signal, respectively, in the cost function, α is the exponent of the fractional order derivative characterizing the dynamics of the heart rate activity y(t), a(t) and b(t) are weighting coefficients for the heart activity and pacing frequency, y min and y max are the minimum and maximum bounds on heart rate activity y(t), y(t i ) is the initial condition, y(t f ) is the final condition, f min and f max are the minimum and maximum allowed bounds on pacing frequency f(t); and
∑
i
=
0
k
(
-
1
)
i
Δ
t
α
(
α
i
)
y
(
(
k
-
i
)
Δ
t
)
-
a
(
k
Δ
t
)
y
(
k
Δ
t
)
+
b
(
k
Δ
t
)
2
z
(
k
Δ
t
)
λ
(
k
Δ
t
)
+
b
(
k
Δ
t
)
β
1
(
k
Δ
t
)
z
(
k
Δ
t
)
-
b
(
k
Δ
t
)
β
2
(
k
Δ
t
)
z
(
k
Δ
t
)
=
0
,
k
=
1
,
…
,
N
∑
i
=
0
N
-
k
(
α
i
)
(
-
1
)
i
λ
(
(
k
+
i
)
Δ
t
)
Δ
t
α
-
w
(
k
Δ
t
)
(
k
Δ
t
)
2
[
y
(
k
Δ
t
)
-
y
ref
(
k
Δ
t
)
]
-
a
(
k
Δ
t
)
λ
(
k
Δ
t
)
+
λ
(
Δ
t
N
)
(
t
f
-
t
i
-
k
Δ
t
)
-
α
Γ
(
1
-
α
)
=
0
,
k
=
N
-
1
,
…
,
0
β
2
(
k
Δ
t
)
-
β
1
(
k
Δ
t
)
-
b
(
k
Δ
t
)
λ
(
k
Δ
t
)
z
(
k
Δ
t
)
-
ξ
1
(
k
Δ
t
)
=
f
min
,
k
=
1
,
…
,
N
β
2
(
k
Δ
t
)
-
β
1
(
k
Δ
t
)
-
b
(
k
Δ
t
)
λ
(
k
Δ
t
)
z
(
k
Δ
t
)
+
ξ
2
(
k
Δ
t
)
=
f
max
,
k
=
1
,
…
,
N
which computes the y and λ values for a predefined set of discretization steps.