Nonlinear reactiondiffusion models of self-organization and deterministic chaos: Theory and possible applications to description of electrical cardiac activity and cardiovascular circulation, Discrete

The paper shows that analytical dynamic models coupled with the available
signal processing methods could be used for describing the
self-organization and chaos degree in the heartbeats propagation and
pressure pulses in ventricular at ejection phase. We proposed a unit
analytical approach that could be associated with real ECG and pressure
pulses signal processing. Our findings confirm that the real-time computer
monitoring of the main cardiovascular parameters obtained by the use of
analytical models and verified by signal processing of real clinical data
may be considered as available method for measuring and controlling
self-organization and chaos degree in pulse propagation.


Introduction
Our main goal is application of nonlinear reaction-diffusion models coupled with realtime ECG and pressure pulses signal processing to monitoring and analysis of the main cardiovascular parameters.For effective applications of the new methods of discrete nonlinear dynamics the continuous models should be adopted by novel approaches to discretizations (see, e.g., Gontar [5]).In the present work the coupled analytical and signal processing approach to this problem is proposed.In the recent years a number of experimental and signal-processing investigations of cardiovascular disorders: atrial and ventricular tachycardia, ventricular fibrillation, acute hypertension, and so forth have been conducted (see, e.g., [3,9,19]).On the other hand, the new concepts of nonlinear dynamics like coherent structures and solitons, fractals and deterministic chaos have been developed (for references see [17]).From our point of view, the above-mentioned concepts coupled with signal processing of biomedical data may be useful in problems of measuring and controlling of self-organization and chaos degree in heartbeat and pulse propagation.We consider the possibility of construction of the so-called hybrid models, the pressure-flow relationship may be approximated by a quadratic function the pressure pulses have been presented as sequences of compact solitary waves (compactons).An analysis of the recurrence energy evolution allows us to conclude that the pulse regularity ranges and critical conditions of transition to full deterministic chaos are mainly determined by the ratio of the elasticity coefficients of the ventricle and aorta walls.This conclusion is in accordance with the investigations of the end-systolic pressure-volume relation in left ventricle (see the main results and references in [16]).The suggested model allows to estimate some important hemodynamic parameters that are not measured directly by conventional methods (e.g., the pulse velocity in the root of the aorta).The obtained quantitative estimation of the critical values of these structural parameters may be applied to the problem of prediction and early diagnosis of transition from arrhytmia to ventricular fibrillation.

Strongly nonlinear model of electrical cardiac activity.
For the simulation of electrical cardiac activity the multidimensional FitzHugh-Nagumo (FHN) system has been considered (see [2]): ( The variable of state u (activator) represents the properties of the membrane potential and excitability, v (inhibitor) is responsible for accommodation and refractoriness.Both α and b are time constants.The parameter c is defined as an adaptive function of the u back time variation δu = u t − u t−∇t .It may be expressed in sigmoidal form: where c 0 and c 1 are constants that scale the response of ∂u/∂t, ∂v/∂t, and δu.During the depolarization ∂u/∂t > 0 → δu > 0 and during repolarization ∂u/∂t < 0 → δu < 0. Although c is a variable parameter, it is kept bounded within the range of 0.2c 0 and 2.2c 0 .
Below is suggested a generalized FitzHugh-Nagumo model in which the diffusion tensor depends on the potential gradient: (2.3) Here the coefficients d i (∇u) of the diffusion are considered as power functions of the heart potential gradient components.From our point of view, it may be associated with self-excitation properties of the cardiac cells depending on electrical potential variations.Below was considered the case of linear inhibitor variation for every heartbeat under the assumption that ∂v/∂t = A n = const.It is supposed that the cardiac cell electrical potential dynamics is governed by modified FHN equation for every heartbeat with coefficients changing from beat to beat.Thus the heartbeats dynamics may be described in the average by a sequence of equation with cubic nonlinearity: , it is reasonable to consider the magnitude β n = A n c − α as a parameter that implements the feedback control of the heart beat-to-beat interaction of activator and inhibitor intensity.We consider solutions of traveling-wave type: (2.5) Here constant V n denotes the pulse velocity during the nth heart beating.Substituting the last expression into (2.4) gives an ordinary differential equation: We denote by D(g) = 3 i=1 g 0 d i (z)dz the function that describes the diffusion potential, D(0) = 0.The structure of this function should be approximated by analysis of experimental electrophysiological data.Evidently, (2.4) may be rewritten as follows: The multiplication of (2.6) by term f n and subsequent integration over the interval [s n, s] (s n ≤ s ≤ s n+1 ) (where s n ,s n+1 are consequence of maximal (peak) points of ECG) gives the following relationship that may be considered as equality that determines the heartbeats "electrical energy" balance: (2.8) The function E n (s) may be considered as the "averaged electrical energy" variability on the interval [s n ,s].This equality coupled with real ECG may be used for the dynamic estimation of the averaged pulse-wave velocity.This may be confirmed by the following motivations.Since f i = 0 for "peak points," the last equation on s = s n+1 for an arbitrary heart beating cycle may be presented as follows: (2.9) V. Kardashov et al. 5 Here A i are maximal amplitudes on ECG (R-values), b 1 = b − 1 and U( f ) may be considered as a potential function of heart electrical activity.The last equality implies an important formula for the averaged heart pulse velocity: It is clear from this formula that the pulse-wave velocity V n is a ratio of the heart potential variation on the R-points of ECG to the mean value of its "kinetic" energy.Since the last formula contains the main structural parameters, V n can be considered as the control parameter that regulates the stability of heartbeats and may cause its instability on an abnormal variation.One can apply the equality (2.9) coupled to real ECG signal for estimation of chaos (or turbulence) degree in heartbeats by direct calculation and analysis of the sequence (E n,n+1 ) N n=1 regularity properties.For such analysis the modern methods of nonlinear discrete dynamics and applied deterministic chaos can be used.Below we present the results of our computer simulation using a real ECG signal.On the other hand, formula (2.10) can be applied to the analysis of the dynamics of the averaged heart potential velocity corresponding to other important parts of ECG, for example, to ST complex.It is known that abnormal alterations of the ST complex are associated with disorders in repolarization process and may be considered as a predictor of transition to ventricular fibrillation (see, e.g., [1,7]).In the presented figures one can observe the fluctuations in the dynamics of the parameter V n and other associated parameters corresponding to the presence of the ventricular tachycardia features in the real ECG.
The above-mentioned relationships determine the dependence of the value V n on the main structural parameters.Of a particular importance is the above-mentioned parameter β n = u n − v n , which in our model is a constant value, changing from one heartbeat to another.It allows us to assume that the pulse-wave velocity (i.e., an important physiological parameter) is controlled, particularly, by the difference of heart electrical potential and inhibitor factors that may be associated with the sympathetic and parasympathetic nervous systems interaction.One can also suppose that this process is implemented by the feedback interaction of the parameters c n and V n .From our point of view, it may help to reveal the nature and physiological mechanisms of the heart rate variability.
From the above discussion follows the importance of constructing the beat-to-beat time series of the above-mentioned general electrophysiological parameters using ECG data.By integration of the equality (2.6) over the interval [s n, s n+1 ] one obtains the relationship: where We consider a diffusion function of the type: 6 Nonlinear reaction-diffusion models Here the exponent m characterizes the nonlinear diffusion intensity.It is evident that when m = 0 the diffusion coefficient is constant.Multiplication of (2.6) in turn by f n , f 2 n and subsequent integration over [s n, s n+1 ] (taking into account the relationships with respect to variables W n , K n , L n , and D n , where The magnitude of c n may be considered as control parameter, responding on the neural feedback regulation by using sympathetic and parasympathetic interaction (in agreement with the variations of the parameter β n ).Thus relationships (2.10)-(2.14)make up the system of linear algebraic equations that allows to construct the time series of the abovementioned cardio-physiological parameters.
For an averaged description of the heart potential diffusion dynamics we propose to calculate the time series of the parameter: On the other hand, our studies implies that if the conditions are fulfilled, which corresponds to the assumption that m > 0, the finite-localized (compact) waves are formed [10].From our point of view, the strong nonlinear phenomena of the steady state localized structure forming in the heart electrical potential distribution can be considered as a powerful interblocking factor that may cause acute heart failure, arrhythmia, and a transition to venticular fibrillation.On the other hand, in recent investigations, the formation of steady state waves in the heart potential propagation is considered as possible, so-called gentle, defibrillation factor (see, e.g., [6]).
In the following figures and tables, one may observe the interactive response of the above-mentioned system of parameters to the presence of the features of the arrhythmia and heart failure in the patients of ECG (see Figure 3.1 and Tables 4.1 and 4.2).It is evident that on the time intervals of the "Holter" signal marked as vt1-vt4 corresponding to presence of ventricular tachycardia in patient heartbeats is observed the reducing of the main nonlinear dynamic and statistical characteristics.

Strongly nonlinear model of pulse propagation in ventricle-aorta system
A multidimensional generalization of the Korteveg-de Vries equation, is considered as an averaged dynamic model of the blood pressure dynamics in the system ventricle-aorta during the ejection phase.Here the function f (u) describes blood flux through the ventricle and g(u) characterizes the blood injection to the aorta, the function u(t,x) describes the output blood pressure distribution.The structural functions f (u),g(u) should be determined by the method of dynamic identification and improved by using experimentally recording signals.
We will consider that pressure-flow relationship (compliance) in ventricle can be described by function of the type: and the averaged injection (blood diffusion) can be described by the function g(u) = du 2 .
Here one may assume that coefficients a, a 1 depend on the geometrical properties of the ventricle and coefficient b depends on its wall elasticity properties.Correspondingly, coefficient d depends on the aorta elasticity (see, e.g., [2]).We consider solutions of (3.1) of the traveling-wave type: Here v = const is the velocity of the traveling wave.By substitution of expression (3.3) into (3.1)one can obtain an ordinary differential equation with respect to the function ϕ(s) (here s = 3 i=1 x i + vt) is a self-similar variable: Integration reduces (3.4) to a second-order equation with an arbitrary number parameter C 0 that is determined by the initial conditions: Assuming that and multiplying (3.5) by expression (d/ds)g(ϕ) = g (ϕ)(dϕ/ds) one may obtain the firstorder equation that determines the total energy conservation law: V. Kardashov et al. 9 Here the function Ψ(ϕ,v) is determined by formula: where the constants E,ϕ 0 denote the total energy and initial value.Furthermore we assume that g(ϕ 0 ) = g (ϕ ) = 0. On the other hand, it is known that if the function g(u) has singular points, satisfying to condition g (u) = 0, then, as evident from equality (3.7), the corresponding solutions are nonanalytic (see [15,11]).These singular solutions can be used to model inelastic shear stress in large arteries.
Formulas (3.7), (3.8) determine the function that gives the dependance of the local momentary energy from structural parameters and the pulse velocity v. Our main assumption is that analytical and computer analysis of this function using the concepts of the modern chaotic vibrations and deterministic chaos theory (see, e.g., [14]) allows a quantitative estimation of the nonregularity and chaos degree of the pressure pulses and allows to find the optimal ranges of variation of the above-mentioned parameters.On the other hand, by the direct substitution into (3.7) the corresponding parameters of the real pressure pulse signal one may construct beat-to-beat time series that describe the dynamics of the associated hemodynamic parameters.It may be applied to computer monitoring and implementation of diagnostic procedures using the modern concepts of nonlinear dynamics (fractal dimension, Lyapunov exponents, etc.).
Below we suppose that ϕ 0 = 0 and s(0,v) = 0.By separating the variables in (3.7) one may obtain the expression that determines the inversion to solution It can be shown that for wide classes of structural functions f ,g, for instance, for even functions, the solutions of problem (3.1) determined by formula (3.9) are by three parameters dependent periodic function.On some critical values of these parameters and structural functions coefficients these solutions transfer to transitional waves (transients) of the two following types: compact solutions of the solitons type (finite-localized transients: compactons or peakons [11,15]) and compact shock waves that, in contrast to the traditional transients, have finite-transitional time or periodic structure (see the recent author's paper [10]).
We suppose that f (u) = a + a 1 u + bu 2 , g(u) = du 2 .Such approximation of the structural functions conforms with the averaged physical conditions of the flow propagation in elastic tube (laminar Poiseuille's flow profiles and experimental displaying form of the pressure-flow relationship).By using the general formula (3.8) one may obtain (3.10) 10 Nonlinear reaction-diffusion models Below we will consider the localized solutions of (3.7) (see the conditions of its formation, e.g., in [10]).It is evident from (3.7) that the conditions of the solution localization are the equality E = 0 and the requirement that the number 0 = 0 is two-multiple root of the equation The nonzero roots of this equation are determined by the formulas Under these conditions there exist two harmonic solutions that according to (3.9) are determined by the formulas (3.13) These harmonic solutions are suggested as a basis for the approximation of pressure and flow pulses in the root of aorta.
The experimentally constructed localized waves may be observed upon restriction of output flow (e.g., at moment of the blood pressure measuring).From the physical point of view, one may assume that this condition can be modeled by offsetting the output flow by initial arbitrary constant, for example, one may assume that It is evident that under these conditions A + dC 0 = 0 the value ϕ = 0 is two-multiple root of the equation Taking this into account, from formulas (3.13) one can obtain the exact expressions for finite localized transients of positive and negative polarities, respectively: It is known from soliton properties that the solitary wave ϕ + (s,v) (of positive polarity) runs from left to right (from heart to peripheral regions) and ϕ − (s,v)-in the opposite direction.This corresponds to the known fact of the existence of the forward and backward pulse waves in the cardiovascular system [18].
Using exact expressions (3.16) for the approximation of experimentally recording signals one can estimate the main structural parameters a i ,b i ,d i and pulse velocity v in the root of aorta.It is reasonable, from our point of view, to approximate the systolic wave of the pulse pressure signal by function ϕ + (s,v) and the dicrotic wave (corresponding to diastolic part of the pulse pressure) by function ϕ − (s,v).Below we give the results of computer simulation of an approximation of the real pulse pressure signal, recorded from the radial artery of the patient by using a piezoelectrical sensor.In Figure 3.2 one can see the results of the root-mean square approximation of one cycle of this pulse pressure signal.
For identification of the main structural parameters the following equalities were used: Here A i , L i denote the nondimensional values of the amplitude and width of the systolic and dicrotic waves that may be measured on the recorded pressure pulse signals, and v i ,b i ,d i are corresponding values of the structural parameters that should be determined.
In addition one may use the relationships where t 1 , t 2 are the time values, corresponding to the maximal amplitudes of the systolic and dicrotic waves.The system of equalities (3.16)- (3.18) allows one to obtain an exact formula for calculating the main structural parameters mentioned above.The most important are the velocities of the dyastolic wave v 1 , dicrotic wave v 2 , and the resulting pressure pulse velocity v: Here the value z corresponds to the distance of the pressure wave from the root of aorta.
The obtained values of these parameters are sufficiently realistic: for the values z = 0.005m and z = 0.013m correspondingly.Here the values t 1 = 12 * 10 −3 s, t 2 = 32 * 10 −3 s, L 1 = 56.4,and L 2 = 137.4were determined by the signal processing of the real pressure pulse signal recorded from the radial artery.It is evident that formula (3.18) is applicable for sufficiently small values of the variable z.From our point of view, at a distance from root of the aorta one should take into account the regulating influence of the peripheral vessels impedance (resistance).For more exact determining of these parameters the correlation coefficients should be introduced.Recent biomedical and clinical investigations conform the crucial importance of these parameters [13].
Below we use a simple "energetic approach" for determining critical relationships on the transition to chaotic dynamics.Since the obtained solitary waves present the so-called separatrix solutions, they determine the bounds of a set of regular (periodic) solutions.They determine also the size of the so-called separatrix layers on the phase plane, generated the chaotic trajectories.
It is important that according to the above-mentioned "energetic" interpretation of formula (3.7) one can see that on the quadratic approximating of the structural functions the momentary kinetic energy in the root of aorta may be presented as the function (3.21) Here λ = L 2 .If we consider the pulse propagation in the vicinity of the aorta as recurrent process of exchanging the pulse pressure and blood flow kinetic energy, the following conclusion may be done: pulse propagation chaos degree may be described by Feigenboum's theory (see [4]).
According to Feigenboum's theory, under the condition 3 < λ < 1 + √ 6 the loss of stability of the pulse energy and the doubling of frequency is displayed.Furthermore, after multiple doubling of frequency under the critical condition pulse transition to complete chaos dynamics may be observed.Thus for the momentary kinetic energy of the compact solutions (3.16) after a sequence of double bifurcations displays the transition to complete chaos (Figure 3.3).On the other hand, according to formulas (3.16) the number L ∞ presents the critical value of the pressure pulse frequency on which the transition to deterministic chaos can be observed.Taking into account these results, one may consider the value L ∞ as a prognostic index, characterizing the pulse transition to deterministic chaos.

Possible applications to biomedicine and conclusions
The main points of our approach are as follows.
(1) Pressure pulse signal, recorded by the noninvasive technique, should be used for the analysis of the main cardiovascular parameters by signal processing procedures.The obtained data may be used for the improvement of the above-mentioned mathematical models.For this purpose the dynamic recognition of the model structural parameters can be used.It may be implemented particularly by application of the experimentally recorded signals and subsequent medical verification.By sufficient number iterations of this process relatively accurate values of the model structural parameters can be obtained.(2) Constructed analytical formulas for the main cardiovascular parameters and exact analytical expressions for solutions, coupled with recorded signals, allow to estimate the above-mentioned hemodynamic and electrophysiological parameters and to implement the real-time monitoring.Obtained time series display the interactive variability of the main cardiovascular parameters that allow to implement the complex analysis.Thus, this approach allows to use for diagnostics of the cardiovascular diseases the exact solutions and analytical formulas coupled with signal processing of the pressure pulse signal.
(3) The method can be applied for estimating the critical bounds of the pressure pulse regularity and chaos degree by using discrete dynamic algorithms.By application of the modern nonlinear dynamic concepts and calculating the main statistical indices (correlation dimension, correlation entropy, fractal dimension, etc.) one can implement the computer assisted diagnostic and prognostic investigations of the cardiovascular system conditions of the patient.Table 4.1 presented the results of pulse-to-pulse monitoring of the averaged dynamics of the above-mentioned index L for number of the pressure pulse signals, recorded by piezoelectrical sensor from radial artery of the patients with different degree of arrhythmia.Table 4.2 presents the corresponding values of the largest lyapunov's exponents (LLE) (see, e.g., [14]).It is known that the range of positive value of the LLE characterizes the chaos degree of the signal.One may see the maximal value of the LLE for the signal, corresponding to File 15.As have been validated by medical diagnosis the corresponding patient condition was classified as pulse arrhythmia.

Figure 3 . 1 .
Figure 3.1.Interactive sensitivity of the main electrophysiological parameters to presence of VT.

Table 4 .
1. Reduction of the main nonlinear dynamic characteristics on the presence of arrhythmia.

Table 4 .
2. Estimation of pressure pulses chaos degree by largest Lyapunov's exponent.