A Novel Nonlinear System Identification for Cerebral Autoregulation in Human: Computer Simulation and Validation

Abstract

Cerebral autoregulation in healthy humans was studied using a novel methodology adapted from Bendat nonlinear analysis technique. A computer simulation of a high-pass filter in parallel with a cubic nonlinearity followed by a low-pass filter was analyzed. A linear system transfer function analysis showed an incorrect estimate of the gain, cut-off frequency, and phase of the high-pass filter. By contrast, using our nonlinear systems identification, yielded the correct gain, cut-off frequency, and phase of the linear system, and accurately quantified the nonlinear system and following low-pass filter. Adding the nonlinear and linear coherence function indicated a complete description of the system. Cerebral blood flow velocity and arterial pressure were measured in six data sets. Application of the linear and nonlinear systems identification techniques to the data showed a high-pass filter, like the linear transfer function, but the gain was smaller. The phase was similar between the two techniques. The linear coherence was low for frequencies below 0.1 Hz but improved by including a nonlinear term. The linear + nonlinear coherence was approximately 0.9 across the frequency bandwidth, indicating an improved description over the linear system analysis of the cerebral autoregulation system.

Appendix

The goal of this appendix is to provide the equations used in this manuscript. We will use G(f) to represent the one-sided frequency response functions computed from standard procedures using fast Fourier transforms of input and output time-domain signals.

Let x(t) represent the input signal, and y(t) the output signal, as illustrated in Fig. 1. As shown, the output signal, y(t) consists of contributions from a linear transformation of the input signal, by the gain and phase of H(f), and a nonlinear transformation modified by the gain and phase of A_n(f). An important consideration is a correlation, or shared energy, between the input signal, x(t), and the input signal as modified by a nonlinearity. In the figure this is indicated by v(t) but here we define it as x₂(t). As a descriptive example of the possible correlation between the input signal and a nonlinearity, let

$$x\left( t \right) = kcos\left( {\omega t} \right),$$

be the input signal where k is the amplitude, and ω = 2πf, where f is frequency. Now if

$$y\left( t \right) = x + x^{3}$$

then

$$y\left( t \right) = kcos\left( {\omega t} \right) + k^{3} cos^{3} \left( {\omega t} \right).$$

Expanding the second term yields

$$k^{3} \left( {\frac{1}{4}\cos \left( {3\omega t} \right) + \frac{3}{4}\cos \left( {\omega t} \right)} \right).$$

Here, in addition to creating the third harmonic, 3ω, the nonlinearity creates energy at the same frequency, ω, as the fundamental frequency. Thus y(t) contains energy at the fundamental frequency from the linear term and the nonlinear term. The goal is to parse the contribution from the two sources.

For random data input, eliminating the contribution of the nonlinearity to the energy of the linear path can be completed in the following manner using conditioned records.3 Let x₁(t) be the signal that will follow the linear path and x₂(t) the signal after transformed by the nonlinearity (i.e. v(t), in Fig. 1) following the nonlinear path. Let G₁₁(f), and G₂₂(f) be the two input auto spectral density functions, and G_yy(f) the output autospectral density function. Here the subscript 11 refers to x₁(t), 22 refers to x₂(t), and yy refers to y(t). G_1y(f) is defined as the cross-spectral density function between x₁(t) and y(t), G_2y(f) is between x₂(t) and y(t), and G₁₂(f) is the cross-spectral density function between x₁(t) and x₂(t). For simplicity of presentation, f is eliminated in the following equations. Define the following frequency response function between the linear and nonlinear input signals, x₁(t) and x₂(t) as,

$$L_{12} = \frac{{G_{12} }}{{G_{11} }}.$$

Remove the influence of the linear path on the spectral density of x₂(t) by

$$G_{22.1} = G_{22} - L_{12} *G_{21}$$

where 22.1 refers to the autospectral density function of x₂(t) with the effects of x₁(t) removed, and G₂₁ is the conjugate of G₁₂.

Compute the optimum frequency response function between the x₁(t) and y(t) as

$$H_{o} = \frac{{G_{1y} }}{{G_{11} }} = L_{1y} .$$

Now the spectral density function between x₂(t) and y(t) with the effects of x₁(t) removed is given as

$$G_{2y.1} = G_{2y} - L_{1y} *G_{21}$$

and the frequency response function for the nonlinear path is

$$L_{2y} = \frac{{G_{2y.1} }}{{G_{22.1} }} .$$

The A systems presented in Fig. 7, are obtained in the following manner,

$$A_{1} = L_{1y} - \left( {L_{12} *L_{2y} } \right)$$

$$A_{2} = L_{2y} .$$

Notice that if there is no correlation between x₁(t) and x₂(t), then L₁₂ = 0 and A₁ = H_o.

The coherence functions representing the proportion of total energy due to the linear and nonlinear path are obtained as follows,

$$L_{c} = \frac{{/L_{1y} /^{2} *G_{11} }}{{G_{yy} }}$$

$$NL_{c} = \frac{{/L_{2y} /^{2} *G_{22.1} }}{{G_{yy} }}$$

where L_c and NL_c indicate linear and nonlinear coherence respectively.