In this appendix, the necessary theoretical underpinning of the fundamental equation 7, which constitutes the basis of the VFR method, is provided.

### A1: fundamental relation between measured voltages and volume changes

In this subsection, the validity of the assumed linearity in eq.(7) is examined. To this end, we first focuson only the ventricles and assume, for the moment, that the atria do not exist, and examine to what extent the measured changes in voltages, as combined into the vector **v**
_{
meas
}, may be considered to be proportional to the volume change *V*(*τ*)−*V*(0). Furthermore, we concentrate on the stroke volume of each heart beat, and enumerate the heart beats using an index (*n*) (*n*∈{1,2,3,…}), and define {\stackrel{~}{V}}_{\left(n\right)} as the stroke volume of beat number *n*. Let **v**
_{
meas,(n)}be defined as the **v**
_{
meas
}vector at the moment in time of maximum volume change during cardiac cycle number *n* (i.e., the end-systolic moment in cardiac cycle number *n*), and let the index (0) refer to the first heartbeat after the start of the measurements (baseline condition, *n*=0). Evidently, in the course of time after *n*=0, during later heartbeats, the stroke value may differ from the initial stroke value {\stackrel{~}{V}}_{\left(0\right)}

In order to evaluate the limitations to the assumed linearity in eq.(7), it is examined to what extent the following “assumption of linearity” \mathcal{L} is valid:

\text{assumption}\phantom{\rule{2.36043pt}{0ex}}\phantom{\rule{2.36043pt}{0ex}}\mathcal{L}\phantom{\rule{2.36043pt}{0ex}}:\phantom{\rule{2.36043pt}{0ex}}\phantom{\rule{2.36043pt}{0ex}}\text{if}\phantom{\rule{2.36043pt}{0ex}}{\stackrel{~}{V}}_{\left(n\right)}\phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}{\lambda}_{\left(n\right)}\phantom{\rule{0.3em}{0ex}}{\stackrel{~}{V}}_{\left(0\right)}\phantom{\rule{2.36043pt}{0ex}}\phantom{\rule{2.36043pt}{0ex}}\text{then}\phantom{\rule{2.36043pt}{0ex}}\phantom{\rule{2.36043pt}{0ex}}{\mathbf{v}}_{\mathrm{meas},\left(n\right)}\phantom{\rule{0.3em}{0ex}}\approx \phantom{\rule{0.3em}{0ex}}{\lambda}_{\left(n\right)}\phantom{\rule{0.3em}{0ex}}{\mathbf{v}}_{\mathrm{meas},\left(0\right)}

(16)

in which *λ*
_{(n)}is a scalar (representing the proportionality factor for cardiac cycle number *n*), and in which the 6 values {\phi}_{\mathrm{meas}}^{\left(k\right)} (*k*∈{1,2,…,6}) in **v**
_{
meas
}are caused by volume changes in the only volumes present, i.e. the ventricles. The assumption \mathcal{L}thus reads that if during cardiac cycle # (n) the stroke volume {\stackrel{~}{V}}_{\left(n\right)} equals *λ*
_{(n)} times the initial stroke volume {\stackrel{~}{V}}_{\left(0\right)}, then all 6 voltages in **v**
_{
meas,(n)}would equal *λ*
_{(n)}times the initial voltages during cardiac cycle # (0) as well.

Evidently, based on the general fact that calculation of voltages from conductivity distributions is a non-linear problem, it is clear on beforehand that the assumption \mathcal{L}will not be true with exact precision; hence the ≈-sign in eq.(16)

In the experiments on human volunteers, an initial calibration procedure using Doppler Ultrasound has been performed at the start of the experiment, and hence the value of the actual initial {\stackrel{~}{V}}_{\left(0\right)} has been provided by Doppler Ultrasound for each volunteer.

The assumption \mathcal{L}can only be valid (approximately) if, for any value of *n*>0, the shape change of the volume, during the volume change during cardiac cycle number *n*, is not completely deviant from the shape change of the volume during *n*=0. In other words, we make the assumption that sudden ruptures or other sudden abnormal patterns of volume changes are excluded. We will refer to this assumption as the “assumption of regularity \mathcal{R}”, and will specify this assumption \mathcal{R} more precisely later on in this subsection.

In order to assess the limits of accuracy of the assumption \mathcal{L}, a fundamental relation between the measured relative voltage changes *φ* and volume changes of blood-filled lumina is derived below, starting with explaining the nature of the applied current field.

The oscillation frequency *f* of the measured AC voltages, which equals the oscillation frequency of the current injected into the patient via the current feeding electrodes on the legs and arms or head, is high in comparison to the physiological processes like respiration or the motion of the heart during the cardiac cycle. Furthermore, the frequency *f* is low enough to let the values of 2*Πfε* (in which *ε* is the dielectric permittivity) be negligible with respect to the differences in specific conductivity *σ* for the relevant tissues [22]. Therefore only the amplitude of the measured AC voltages needs to be considered; in our set-up each {\Phi}_{\mathrm{meas}}^{\left(k\right)}\left(\tau \right) value was defined as the *amplitude of the oscillating voltage* measured between the electrodes *e*
_{
k
} and *e*
_{
k−1} at time *τ*. This amplitude of the voltage was established every 5 milliseconds.

Directly below, a relation between the measured voltage changes *φ* and volume changes of blood-filled lumina is derived, on basis of the law of conservation of charge (the continuity equation) and Ohm’s law.

Since the current sources are applied to the head or arms and legs, there are no current injection electrodes present on the thorax. Therefore, the continuity equation for points inside the thorax reads:

\nabla \xb7\mathbf{J}\left(\mathbf{x}\right)=\phantom{\rule{1em}{0ex}}0\phantom{\rule{1em}{0ex}}\text{for all}\phantom{\rule{1em}{0ex}}\mathbf{x}\in \mathcal{T}

(17)

in which **x** is a three-dimensional spatial coordinate ( i.e. **x** denotes a point within 3D space), and **J**(**x**) denotes the current density at point **x**, and \mathcal{T} denotes the complete thorax, including the thoracic skin, and everything inside the thorax.

Substituting Ohms law **J**=*σ* **E** in eq. (17), in which *σ* denotes the conductivity and **E** denotes the electric field, and applying the general mathematical identity

\nabla \xb7\left(p\mathbf{Q}\right)=p\nabla \xb7\mathbf{Q}+\mathbf{Q}\xb7\nabla \phantom{\rule{1em}{0ex}}p\phantom{\rule{1em}{0ex}}\text{for any scalar}p\text{and vector}\phantom{\rule{1em}{0ex}}\mathbf{Q}

we may rewrite eq. (17) as

\sigma \nabla \xb7\mathbf{E}+\phantom{\rule{2.36043pt}{0ex}}\mathbf{E}\xb7\nabla \sigma =0

(18)

Using \frac{1}{\sigma}\nabla \sigma =\nabla log\sigma, this can be written as

\nabla \xb7\mathbf{E}=-\phantom{\rule{0.3em}{0ex}}\mathbf{E}\xb7\nabla log\sigma

(19)

Let the 3D vector field −∇log(*σ*(**x**))be denoted as **s**(**x**).

Therefore we now have:

\nabla \xb7\mathbf{E}\left(\mathbf{x}\right)=\mathbf{E}\left(\mathbf{x}\right)\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\mathbf{s}\left(\mathbf{x}\right)

(20)

This is an implicit equation for **E**(**x**) for a given **s**(**x**) distribution. For a model of the human anatomy in which a limited number of tissue types are defined (such as “blood”, “muscle tissue”, “bone”, etc.), the spatial conductivity distribution *σ*(**x**) is piecewise constant; i.e. the value of **s**(**x**)is zero everywhere, except at the boundary between two different tissue types. For a sharp boundary, we have [21]: \mathbf{s}\left(\mathbf{x}\right)=2\widehat{\mathbf{n}}\frac{{\sigma}_{1}-{\sigma}_{2}}{{\sigma}_{1}+{\sigma}_{2}}, in which \widehat{\mathbf{n}} is the normal unit vector perpendicular to the boundary surface at point **x**, and *σ*
_{1}and *σ*
_{2} are the specific conductivities of the two tissue types that meet at the boundary surface.

Since \nabla \xb7\mathbf{E}=\frac{\rho}{\epsilon} (according to the Maxwell equations, in which *ρ* is the charge density and *ε* is the dielectric constant), the right-hand side of eq.(20), i.e. **E**(**x**) · **s**(**x**), may be regarded as a charge density divided by the dielectric permittivity *ε*. Therefore, points within the thorax having non-zero values of **s**(**x**)may be considered as equivalent “secondary” sources, responsible for non-zero divergence of **E** within the thorax at points on boundary surfaces. The effect of these secondary sources on the rest of the thorax (including the measuring electrodes on the thoracic skin) may be calculated using the Green function *G* of the Poisson equation.

In order to use this for the particular case of ventricular volume changes, consider eq.(20) at two instances in time:∇·**E**
_{0}(**x**)=**E**
_{0}(**x**) · **s**
_{0}(**x**) at time *τ*=0, i.e. at the time of the R-peak of the ECG, and ∇·**E**(**x**,*τ*)=**E**(**x**,*τ*) · **s**(**x**,*τ*)at some later instant during the same cardiac cycle (*τ*≠0). Subtracting the equation for *τ*=0from the equation for *τ*≠0, and using the extra symbols \stackrel{~}{\mathbf{E}} and \stackrel{~}{\mathbf{s}} with \stackrel{~}{\mathbf{E}}(\mathbf{x},\tau )=\mathbf{E}(\mathbf{x},\tau )-{\mathbf{E}}_{0}\left(\mathbf{x}\right) and \stackrel{~}{\mathbf{s}}(\mathbf{x},\tau )=\mathbf{s}(\mathbf{x},\tau )-{\mathbf{s}}_{0}\left(\mathbf{x}\right) to indicate the difference fields, yields

\nabla \xb7\stackrel{~}{\mathbf{E}}(\mathbf{x},\tau )=\underset{(\star )}{\underset{\u23df}{{\mathbf{E}}_{0}\left(\mathbf{x}\right)\xb7\stackrel{~}{\mathbf{s}}(\mathbf{x},\tau )}}+\underset{(\star \star )}{\underset{\u23df}{\stackrel{~}{\mathbf{E}}(\mathbf{x},\tau )\xb7\mathbf{s}(\mathbf{x},\tau )}}

(21)

This is still an implicit equation from which \stackrel{\mathbf{~}}{\mathbf{E}} can not be calculated directly.

In eq.(21), the first term ({\mathbf{E}}_{0}\left(\mathbf{x}\right)\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\stackrel{~}{\mathbf{s}}(\mathbf{x},\tau )), indicated by (⋆), represents the first order approximation, whereas \stackrel{~}{\mathbf{E}}(\mathbf{x},\tau )\xb7\mathbf{s}(\mathbf{x},\tau ), indicated by (⋆⋆), represents the self-reference in this equation, which leads to higher-order corrections in the form of a perturbation series if an exact expression for \stackrel{\mathbf{~}}{\mathbf{E}} is derived. In the appendix, the exact solution to eq.(21) is presented in the form of such a perturbation series.

Since the purpose of this subsection is to determine to what degree of accuracy the actual volume change *V*(*τ*)−*V*(0) may be considered to be proportional to the *ψ*(*τ*)(as calculated from the measurements using eq. (11)), we first examine this proportionality on basis of the first term ({\mathbf{E}}_{0}\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\stackrel{~}{\mathbf{s}}\left(\tau \right)) alone, and, subsequently, calculate the higher-order deviations from this proportionality on basis of the perturbation series that constitutes the solution to eq.(21) as a whole.

Using the first order approximation of eq.(21), i.e. \nabla \xb7\stackrel{~}{\mathbf{E}}\left(\tau \right)={\mathbf{E}}_{0}\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\stackrel{~}{\mathbf{s}}\left(\tau \right), the expression for \stackrel{\mathbf{~}}{\mathbf{E}}(\mathbf{x},\tau ) reads:

\stackrel{\mathbf{~}}{\mathbf{E}}(\mathbf{x},\tau )=\nabla [G\phantom{\rule{0.3em}{0ex}}\ast \phantom{\rule{0.3em}{0ex}}\left(\stackrel{~}{\mathbf{s}}(\mathit{\xi},\tau )\xb7{\mathbf{E}}_{0}\left(\mathit{\xi}\right)\right)]

(22)

in which ∗ denotes the 3D spatial convolution, and *ξ* denotes points in 3D space over which the convolution takes place, and *G* denotes the Green function of the Poisson equation (G(\mathbf{x}-\mathit{\xi})=1/\left(4\Pi \sqrt{(\mathbf{x}-\mathit{\xi})\xb7(\mathbf{x}-\mathit{\xi}}\right))). Hence the change in the voltage difference *Φ*^{(k)}between two electrodes *e*
_{
k
}and *e*
_{
k−1} (at positions **x**
_{
k
} and **x**
_{
k−1} on the skin) as function of the time *τ* during one single heartbeat is calculated as:

{\Phi}^{\left(k\right)}\left(\tau \right)-{\Phi}^{\left(k\right)}\left(0\right)={\int}_{\mathrm{thorax}}\phantom{\rule{0.3em}{0ex}}{d}^{3}\xi \left[G({\mathbf{x}}_{k-1}-\mathit{\xi})-G({\mathbf{x}}_{k}-\mathit{\xi})\right]\stackrel{~}{\mathbf{s}}(\mathit{\xi},\tau )\xb7{\mathbf{E}}_{0}\left(\mathit{\xi}\right)

(23)

in which *∫*
_{
thorax
}
*d*^{3}
*ξ* denotes integration over 3D space. with *ξ* denoting points in 3D space. From this equation, a basic relation between voltage changes and volume changes can be derived, as is indicated in Figure 8.

Consider, as a simplified example, the boundary of a single volume of blood within an environment of tissue. This boundary is rendered in Figure 8a for the initial situation at *τ*=0. The value of **s**(**x**,0)is zero everywhere, except at the boundary. At a later instant *τ*
_{
ES
}during the same cardiac cycle (in which the subscript “*ES*” refers to “end-systole”), the volume of blood has reached its maximum volume change.

As is explained in the legend of Figure 8, the subtraction of Figure 8a from Figure 8b yields Figure 8c, in which the equivalent charge distributions \stackrel{~}{\mathbf{s}}\xb7{\mathbf{E}}_{0}\left(\mathbf{x}\right) are rendered that give rise to *Φ*^{(k)}(*τ*
_{
ES
})−*Φ*^{(k)}(0) in eq.(23). All situations rendered in Figure 8 are from one single heartbeat, i.e. cardiac cycle number *n*. The relation with maximum volume change (i.e., the stroke volume {\stackrel{~}{V}}_{\left(n\right)} of cardiac cycle number *n*) is now indicated in Figure 8d. The volume change {\stackrel{~}{V}}_{\left(n\right)}=V\left({\tau}_{\mathrm{ES}}\right)-{V}_{0}, in which *V*(*τ*
_{
ES
}) relates to Figure 8b, and *V*
_{0}relates to Figure 8a, can be calculated as the sum of the set of small grey rectangular beams in Figure 8d. Let the small grey beams be enumerated using the index *i*. Since the direction of the small grey beams, and hence of *δ* **A**
_{
i
}, may be chosen freely when calculating a volume, we choose *δ* **A**
_{
i
}to be parallel to {\hat{\mathbf{e}}}_{0}, in which {\hat{\mathbf{e}}}_{0} is the direction parallel to **E**
_{0}, with

{\mathbf{E}}_{0}={E}_{0}{\hat{\mathbf{e}}}_{0}\text{and}\delta {\mathbf{A}}_{i}={A}_{i}{\hat{\mathbf{e}}}_{0}\text{and}{\mathbf{h}}_{i}={h}_{i}{\hat{\mathbf{e}}}_{0}

Therefore the volume change now reads

{\stackrel{~}{V}}_{\left(n\right)}=V\left({\tau}_{\mathrm{ES}}\right)-V\left(0\right)=\sum _{i}{\mathbf{h}}_{i}\xb7\left(\delta {\mathbf{A}}_{i}\right)=\sum _{i}{h}_{i}{A}_{i}

(24)

The model depicted in Figure 8 can incorporate volume changes in a direction perpendicular to {\hat{\mathbf{e}}}_{0} as well: If the volume is becoming smaller in a direction perpendicular to {\hat{\mathbf{e}}}_{0} , then some rods near the boundary become gradually shorter until they vanish completely (when the length of the rod becomes zero). If the volume is becoming larger in a direction perpendicular to {\hat{\mathbf{e}}}_{0} , some vanished rods with length zero obtain non-zero length again (effectively adding new rods at the boundary of the volume).

Similarly, the charges **s**·**E**
_{0}(**x**)may be considered to be a set of dipoles {\mathbf{p}}_{i}={q}_{i}{h}_{i}{\hat{\mathbf{e}}}_{0} (as explained in Figure 8), enumerated by the index *i*, with dipole length *h*
_{
i
} and charge *q*
_{
i
}.

We will now examine the proportionality between the volume change {\stackrel{~}{V}}_{\left(n\right)}=V\left({\tau}_{\mathrm{ES}}\right)-V\left(0\right) and the voltage differences measured on the thorax resulting from the set of dipoles {\mathbf{p}}_{i}={q}_{i}{h}_{i}{\hat{\mathbf{e}}}_{0}.

At this point, we provide a more precise definition of the “assumption \mathcal{R}” (the assumption of regularity in the volume changes) that has been mentioned at the beginning of this subsection:

\text{assumption}\mathcal{R}:\mathrm{\text{if}}{\stackrel{~}{V}}_{\left(n\right)}={\lambda}_{\left(n\right)}{\stackrel{~}{V}}_{\left(0\right)},\text{then}\forall i:{h}_{i,\left(n\right)}={\lambda}_{\left(n\right)}{h}_{i,\left(0\right)}

(25)

With this assumption \mathcal{R}, and the dipole modeling described above, the voltage change *Φ*^{k}(*τ*
_{
ES
})−*Φ*^{k}(0)in eq.(23) for cardiac cycle *#n* may be written as:

{\Phi}_{\left(n\right)}^{k}\left({\tau}_{\mathrm{ES}}\right)-{\Phi}_{\left(n\right)}^{k}\left(0\right)=\sum _{i}\left[{\mathbf{G}}_{\mathrm{dipole}}({\mathbf{x}}_{k}-{\mathit{\xi}}_{i})-{\mathbf{G}}_{\mathrm{dipole}}({\mathbf{x}}_{k-1}-{\mathit{\xi}}_{i})\right]\xb7{\hat{\mathbf{e}}}_{0}{h}_{i,\left(n\right)}{q}_{i}

(26)

in which {\mathbf{G}}_{\mathrm{dipole}}\left(\mathbf{x}\right)=\widehat{\mathbf{x}}/\left(4\mathrm{\Pi \epsilon}\right(\mathbf{x}\xb7\mathbf{x}\left)\right).

Substituting *h*
_{
i,(n)}=*λ*
_{(n)}
*h*
_{
i,(0)}from eq. (25) and *q*
_{
i
}=2*ε* *E*
_{0}
*A*
_{
i
}(*σ*
_{
blood
}−*σ*
_{
tissue
})/(*σ*
_{
blood
} + *σ*
_{
tissue
})[21] in the simple example of Figure 8, this yields :

\begin{array}{ll}\phantom{\rule{-20.0pt}{0ex}}{\Phi}_{\left(n\right)}^{k}\left({\tau}_{\mathrm{ES}}\right)-{\Phi}_{\left(n\right)}^{k}\left(0\right)& =\sum _{i}\left[{\mathbf{G}}_{\mathrm{dipole}}\left({\mathbf{x}}_{k}-{\mathit{\xi}}_{i}\right)-{\mathbf{G}}_{\mathrm{dipole}}({\mathbf{x}}_{k-1}-{\mathit{\xi}}_{i})\right]\xb7{\hat{\mathbf{e}}}_{0}{h}_{i,\left(n\right)}\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{1em}{0ex}}\times \left(2{A}_{i}{E}_{0}\frac{{\sigma}_{\mathrm{blood}}-{\sigma}_{\mathrm{tissue}}}{{\sigma}_{\mathrm{blood}}+{\sigma}_{\mathrm{tissue}}}\right)=2{\lambda}_{\left(n\right)}{E}_{0}\frac{{\sigma}_{\mathrm{blood}}-{\sigma}_{\mathrm{tissue}}}{{\sigma}_{\mathrm{blood}}+{\sigma}_{\mathrm{tissue}}}\sum _{i}{h}_{i,\left(0\right)}{A}_{i}\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{1em}{0ex}}\times \left[{\mathbf{G}}_{\mathrm{dipole}}({\mathbf{x}}_{k}-{\mathit{\xi}}_{i})-{\mathbf{G}}_{\mathrm{dipole}}({\mathbf{x}}_{k-1}-{\mathit{\xi}}_{i})\right]\xb7{\hat{\mathbf{e}}}_{0}\phantom{\rule{2em}{0ex}}\end{array}

(27)

from which it becomes apparent that, having used the assumption of regularity \mathcal{R}(eq. (25)) and the first order approximation in eq.(22), {\Phi}_{\left(n\right)}^{k}\left({\tau}_{\mathrm{ES}}\right)-{\Phi}_{\left(n\right)}^{k}\left(0\right)is indeed proportional to {\lambda}_{\left(n\right)}\sum _{i}{h}_{i,\left(0\right)}{A}_{i}, and hence assumption \mathcal{L}holds in this case. Furthermore, eq.(27) shows that {\Phi}_{\left(n\right)}^{k}\left({\tau}_{\mathrm{ES}}\right)-{\Phi}_{\left(n\right)}^{k}\left(0\right)is proportional to the incident field strength *E*
_{0}as well.

As has been indicated before however, it is clear on beforehand that the assumption \mathcal{L} will not be true with exact precision, due to the fact that the relation between volumes and voltages in fundamentally non-linear. The strength of the non-linear effects will be examined in the next subsection (A2).

### A2: non-linear effects

In the subsection above, the first order explicit equation (22) has been derived from the general implicit equation (21), by neglecting the self-referential term indicated by (⋆⋆)in eq.(21). The exact version of equation (22) however, i.e. without neglecting the self-referential term indicated by (⋆⋆) in eq.(21), reads:

\stackrel{\mathbf{~}}{\mathbf{E}}\left(\mathbf{x}\right)=\nabla {\int}_{\mathrm{thorax}}\phantom{\rule{0.3em}{0ex}}{d}^{3}\mathrm{\xi G}(\mathbf{x}-\mathit{\xi})\left({\mathbf{E}}_{0}\left(\mathit{\xi}\right)\xb7\stackrel{~}{\mathbf{s}}\left(\mathit{\xi}\right)+\mathbf{s}\left(\mathit{\xi}\right)\xb7\stackrel{\mathbf{~}}{\mathbf{E}}\left(\mathit{\xi}\right)\right)

(28)

and hence the voltage difference between the two electrodes *e*
_{
k
}and *e*
_{
k−1}equals:

\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{\Phi}^{\left(k\right)}\left(\tau \right)-{\Phi}^{\left(k\right)}\left(0\right)={\int}_{\mathrm{thorax}}\phantom{\rule{0.3em}{0ex}}{d}^{3}\xi \left[G({\mathbf{x}}_{k-1}-\mathit{\xi})-G({\mathbf{x}}_{k}-\mathit{\xi})\right]\left({\mathbf{E}}_{0}\left(\mathit{\xi}\right)\xb7\stackrel{~}{\mathbf{s}}\left(\mathit{\xi}\right)+\mathbf{s}\left(\mathit{\xi}\right)\xb7\stackrel{\mathbf{~}}{\mathbf{E}}\left(\mathit{\xi}\right)\right)

(29)

In order to transform equation (28) into an explicit equation for \stackrel{\mathbf{~}}{\mathbf{E}} (and hence solving eq.(29) as well), an iterative substitution method is used that is comparable to the iterative substitution that produces the Born-Neumann series from the Lippmann-Schwinger equation, i.e.: the symbol \stackrel{\mathbf{~}}{\mathbf{E}} appearing in the right-hand side of eq.(28) is replaced by the expression that is formed by the entire right-hand side of eq.(28). This results in a new equation in which \stackrel{\mathbf{~}}{\mathbf{E}} appears again near the end of the right side of the resulting equation, and hence the substitution procedure can be repeated, which yields an iterative substitution procedure.

In the top line of eq.(30) the result of the iterative substitution is rendered, and subsequently “expanded” by distribution through parentheses (in the second line of eq.(30)), yielding a perturbation series:

\begin{array}{ll}\stackrel{\mathbf{~}}{\mathbf{E}}\phantom{\rule{2.36043pt}{0ex}}=& \phantom{\rule{2.36043pt}{0ex}}\nabla G\phantom{\rule{0.3em}{0ex}}\ast \phantom{\rule{0.3em}{0ex}}\left({\mathbf{E}}_{0}\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\stackrel{~}{\mathbf{s}}\phantom{\rule{2.36043pt}{0ex}}+\phantom{\rule{2.36043pt}{0ex}}\mathbf{s}\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\nabla G\phantom{\rule{0.3em}{0ex}}\ast \phantom{\rule{0.3em}{0ex}}\left({\mathbf{E}}_{0}\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\stackrel{~}{\mathbf{s}}\phantom{\rule{2.36043pt}{0ex}}+\phantom{\rule{2.36043pt}{0ex}}\mathbf{s}\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\nabla G\phantom{\rule{0.3em}{0ex}}\ast \phantom{\rule{0.3em}{0ex}}\left(\dots \dots \right)\right)\right)\phantom{\rule{2em}{0ex}}\\ =& \phantom{\rule{2.36043pt}{0ex}}\nabla G\phantom{\rule{0.3em}{0ex}}\ast \phantom{\rule{0.3em}{0ex}}\phantom{\rule{2.36043pt}{0ex}}\underset{\mathcal{A}}{\underset{\u23df}{\left(\phantom{\rule{2.36043pt}{0ex}}\sum _{\vartheta =0}^{\infty}{\left[\mathbf{s}\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\nabla G\phantom{\rule{0.3em}{0ex}}\ast \phantom{\rule{0.3em}{0ex}}\right]}^{\vartheta}\phantom{\rule{2.36043pt}{0ex}}\right)}}\phantom{\rule{2.36043pt}{0ex}}\phantom{\rule{2.36043pt}{0ex}}({\mathbf{E}}_{0}\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\stackrel{~}{\mathbf{s}})\phantom{\rule{2em}{0ex}}\end{array}

(30)

in which, again, the symbol ∗denotes the convolution operator in three dimensions.

The voltage difference between the two electrodes *e*
_{
k
}and *e*
_{
k−1} hence reads:

\phantom{\rule{-20.0pt}{0ex}}{\Phi}^{\left(k\right)}\left(\tau \right)-{\Phi}^{\left(k\right)}\left(0\right)\phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}{\int}_{\mathrm{thorax}}\phantom{\rule{0.3em}{0ex}}{d}^{3}\xi \phantom{\rule{2.36043pt}{0ex}}\phantom{\rule{2.36043pt}{0ex}}\phantom{\rule{0.3em}{0ex}}\left[G({\mathbf{x}}_{k-1}-\mathit{\xi})-G({\mathbf{x}}_{k}-\mathit{\xi})\right]\phantom{\rule{0.3em}{0ex}}\underset{\mathcal{A}}{\underset{\u23df}{\left(\phantom{\rule{2.36043pt}{0ex}}\sum _{\vartheta =0}^{\infty}{\left[\mathbf{s}\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\nabla G\phantom{\rule{0.3em}{0ex}}\ast \phantom{\rule{0.3em}{0ex}}\right]}^{\vartheta}\phantom{\rule{2.36043pt}{0ex}}\right)}}\phantom{\rule{2.36043pt}{0ex}}\phantom{\rule{2.36043pt}{0ex}}({\mathbf{E}}_{0}\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\stackrel{~}{\mathbf{s}})

(31)

An important parameter for the assessment of convergence of the series \mathcal{A}in eq. (30) and eq.(31)) is the fraction \frac{\mid \stackrel{\mathbf{~}}{\mathbf{E}}\left(\mathbf{x}\right)\mid}{\mid {\mathbf{E}}_{0}\left(\mathbf{x}\right)\mid}.

If the fraction \frac{\mid \stackrel{\mathbf{~}}{\mathbf{E}}\left(\mathbf{x}\right)\mid}{\mid {\mathbf{E}}_{0}\left(\mathbf{x}\right)\mid} is small compared to one, then the higher order (large *ϑ*) terms in (31) become rapidly irrelevant with increasing *ϑ*.

Since the cardiac motion causes only a small distortion of the applied current field, we have that \frac{\mid \stackrel{\mathbf{~}}{\mathbf{E}}\left(\mathbf{x}\right)\mid}{\mid {\mathbf{E}}_{0}\left(\mathbf{x}\right)\mid} is indeed small compared to one, as can be shown for the theoretical case of a single spherical volume representing a ventricle in a homogeneous surroundings. For this theoretical case, the maximum value of \frac{\mid \stackrel{\mathbf{~}}{\mathbf{E}}\left(\mathbf{x}\right)\mid}{\mid {\mathbf{E}}_{0}\left(\mathbf{x}\right)\mid} outside the sphere can be calculated exactly, and equals 0.14 for the conductivity values of blood and muscle tissue, and is reached near the surface of the sphere on the central axis of the sphere parallel to the incident field **E**
_{0}. For most other positions **x**, the value of \frac{\mid \stackrel{\mathbf{~}}{\mathbf{E}}\left(\mathbf{x}\right)\mid}{\mid {\mathbf{E}}_{0}\left(\mathbf{x}\right)\mid} is much lower.

We now use eq.(31) to examine the non-linear character of the relation between volume changes and the measured *Φ*^{(k)}(*τ*) − *Φ*^{(k)}(0) by adding various anatomical and cardiological features, starting with the above mentioned simple theoretical case of a single spherical volume of blood (with *σ*
_{
blood
}) inside an infinite homogenous surroundings (with specific conductivity *σ*
_{
tissue
}). From the literature [21] it is known that, in this specific case, the electric field in the surroundings of the sphere is exactly equal to the sum of the incident field plus the field from a single dipole located at the center of the sphere.

As a result, the *Φ*^{(k)}(*τ*) − *Φ*^{(k)}(0)will be exactly proportional to the volume of the sphere.

Non-linear effects however do occur in \mathcal{A} in eq.(31) if other objects are present besides the single sphere. Evidently, the true deviation from linearity in a complex realistic anatomical model can not be rendered in an explicit analytical expression. However, it is possible to derive analytically a set of boundaries that confine the possible deviations from linearity to a well-defined and narrow range for each specific anatomical and cardiological feature. In order to do so, these specific anatomical and cardiological features are translated into sets of spheres (as will be explained below). First, let the change in the incident field on a sphere *S*, caused by the presence of *another* sphere *S*’, be denoted by \left(\right)close="">\n \n \n \n E\n \n \u0306\n \n (\n x\n )\n \n. In order to calculate the confines of the range of the possible deviations from linearity, the following step is taken: for each point on the surface of *S*, the vector \left(\right)close="">\n \n \n \n E\n \n \u0306\n \n (\n x\n )\n \n of the incident field on *S*, due to the presence of the other sphere *S*’, is replaced by the same vector \left(\right)close="">\n \n \n \n \n \n E\n \n \u0306\n \n \n \n max\n \n \n \n, in which \left(\right)close="">\n \n \n \n \n \n E\n \n \u0306\n \n \n \n max\n \n \n \n is the maximum (“worst case”) value of all original \left(\right)close="">\n \n \n \n E\n \n \u0306\n \n (\n x\n )\n \n vectors on the surface of sphere *S*. This produces a stronger non-linearity than in the case of an exact calculation, and hence provides an analytical expression for the upper limit to the strength of the non-linear effects.

The following three fundamental anatomical and cardiological features are considered that each produce a specific deviation from linearity (see Figure 9):

(a) The non-spherical shape of the ventricles (Figure 9a);

(b) The proximity of the heart to a boundary surface (skin-air), viz. the frontal surface of the thoracic skin (Figure 9b);

(c) The cardiological phenomenon that, e.g. during the rapid filling phase of the ventricles, the volume changes in the atria tend to be comparable to the volume changes in the ventricles, but have the opposite sign. In a completely linear model, neglecting all higher order terms (i.e., neglecting the term indicated by (⋆⋆) in eq.(21) ), the volume changes in the atria would contribute to *only* *α*(*τ*)in eq.(7), and the volume changes in the ventricles would contribute to *only* *ψ*(*τ*). However, indirectly, a volume change in the atria will change also the way that a volume change in the ventricles is perceived by the measuring electrodes, because a volume change in the atria will change the strength of the *incident field* on the ventricles, giving rise to a second order effect. (Figure 9c)

As an example, consider the following calculation of an upper limit to the second order deviation from linearity for feature (b). In Figure 9b, the feature of the proximity of the heart to a boundary surface (skin-air), i.e. feature (b), is rendered schematically. The method of image charges [21] is used to replace the skin-air boundary by a mirror image of the heart in a homogenous continuum.

In the spherical modeling, increasing the stroke volume (as described by \stackrel{~}{\mathbf{s}} ) with a factor *λ* is equivalent to replacing \stackrel{~}{\mathbf{s}} with \lambda \stackrel{~}{\mathbf{s}} in its effect on the measuring electrodes. Furthermore, the \left(\right)close="">\n \n \n \n E\n \n \u0306\n \n (\n x\n )\n \n near sphere *S*, caused by the presence of the mirror sphere *S*’, is proportional to the volume of the mirror sphere *S*’, and therefore also proportional to *λ*, and, as a result, the maximum \left(\right)close="">\n \n \n \n \n \n E\n \n \u0306\n \n \n \n max\n \n \n \n is proportional to *λ* as well. Since \left(\right)close="">\n \n \u2223\n \n \n E\n \n \u0306\n \n (\n x\n )\n \u2223\n \n is also proportional to ∣**E**
_{0}∣, the \left(\right)close="">\n \n \u2223\n \n \n \n \n E\n \n \u0306\n \n \n \n max\n \n \n \u2223\n \n may be written as *λη* *E*
_{0}, in which *η* depends on the distance between the spheres. As a result, on basis of eq.(29), the expression {\left[{\Phi}_{\left(n\right)}^{\left(k\right)}\left({\tau}_{\mathrm{ES}}\right)-{\Phi}_{\left(n\right)}^{\left(k\right)}\left(0\right)\right]}_{max}, representing a “worst case” estimation for \left(\right)close="">\n \n \n \n \Phi \n \n \n (\n n\n )\n \n \n (\n k\n )\n \n \n (\n \n \n \tau \n \n \n ES\n \n \n )\n \u2212\n \n \n \Phi \n \n \n (\n n\n )\n \n \n (\n k\n )\n \n \n (\n 0\n )\n \n that has a maximum deviation from linearity, yields:

\begin{array}{ll}{\left[{\Phi}_{\left(n\right)}^{\left(k\right)}\left({\tau}_{\mathrm{ES}}\right)-{\Phi}_{\left(n\right)}^{\left(k\right)}\left(0\right)\right]}_{max}& ={\int}_{\mathrm{thorax}}{d}^{3}\xi \phantom{\rule{2.36043pt}{0ex}}\phantom{\rule{2.36043pt}{0ex}}\left[G\left({\mathbf{x}}_{k-1}-\mathit{\xi}\right)-G\left({\mathbf{x}}_{k}-\mathit{\xi}\right)\right]\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{1em}{0ex}}\times \left({\mathbf{E}}_{0}\left(\mathit{\xi}\right)\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\left({\lambda}_{\left(n\right)}\stackrel{~}{\mathbf{s}}\right(\mathit{\xi}\left)\right)+\left({\mathbf{s}}_{0}\right(\mathit{\xi})+{\lambda}_{\left(n\right)}\stackrel{~}{\mathbf{s}}(\mathit{\xi}\left)\right)\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}{\stackrel{\u0306}{\mathbf{E}}}_{max}\right.\phantom{\rule{2em}{0ex}}\\ ={\int}_{\mathrm{thorax}}{d}^{3}\xi \left[G({\mathbf{x}}_{k-1}-\mathit{\xi})\phantom{\rule{0.3em}{0ex}}-\phantom{\rule{0.3em}{0ex}}G({\mathbf{x}}_{k}-\mathit{\xi})\right]\left(\underset{\propto \phantom{\rule{0.3em}{0ex}}{\lambda}_{\left(n\right)}{E}_{0}}{\underset{\u23df}{{\mathbf{E}}_{0}\left(\mathit{\xi}\right)\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\left({\lambda}_{\left(n\right)}\stackrel{~}{\mathbf{s}}\right(\mathit{\xi}\left)\right)}}\right.\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{1em}{0ex}}+\left(\right)close=")">\underset{\propto \phantom{\rule{0.3em}{0ex}}\eta {\lambda}_{\left(n\right)}{E}_{0}}{\underset{\u23df}{{\mathbf{s}}_{0}\left(\mathit{\xi}\right)\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}{\stackrel{\u0306}{\mathbf{E}}}_{max}}}+\underset{\propto \phantom{\rule{0.3em}{0ex}}\eta {\lambda}_{\left(n\right)}^{2}{E}_{0}}{\underset{\u23df}{{\lambda}_{\left(n\right)}\stackrel{~}{\mathbf{s}}\left(\mathit{\xi}\right)\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}{\stackrel{\u0306}{\mathbf{E}}}_{max}}}& \phantom{\rule{2em}{0ex}}\end{array}\n

(32)

in which the last term, indicated by the third horizontal brace at the far right, is proportional to the *square* of *λ*
_{(n)}, because *λ*
_{(n)} and \left(\right)close="">\n \n \n \n \n \n E\n \n \u0306\n \n \n \n max\n \n \n \n are both proportional to *λ*
_{(n)}. As a result, the \eta {\lambda}_{\left(n\right)}^{2}{E}_{0} under the third horizontal brace at the far right represents the second-order effect, which is a deviation from the first order linear relation between \left(\right)close="">\n \n \n \n \Phi \n \n \n (\n n\n )\n \n \n (\n k\n )\n \n \n (\n \n \n \tau \n \n \n ES\n \n \n )\n \u2212\n \n \n \Phi \n \n \n (\n n\n )\n \n \n (\n k\n )\n \n \n (\n 0\n )\n \n and *λ*
_{(n)}.

Since {\phi}_{\left(n\right)}^{k}\left({\tau}_{\mathrm{ES}}\right)=\left({\Phi}_{\left(n\right)}^{\left(k\right)}\right({\tau}_{\mathrm{ES}})-{\Phi}_{\left(n\right)}^{\left(k\right)}(0\left)\right)/{\Phi}_{\left(n\right)}^{\left(k\right)}\left(0\right)(by definition, see eq.(2)) and \left({\Phi}_{\left(n\right)}^{\left(k\right)}\right({\tau}_{\mathrm{ES}})-{\Phi}_{\left(n\right)}^{\left(k\right)}(0\left)\right)\ll {\Phi}_{\left(n\right)}^{\left(k\right)}\left(0\right)(empirical fact, see e.g. [15] and [23]), and {\Phi}_{\left(n\right)}^{\left(k\right)}\left(0\right)\propto {E}_{0}, the three proportionalities indicated under the braces in eq.(32), i.e. \propto \phantom{\rule{0.3em}{0ex}}{\lambda}_{\left(n\right)}{E}_{0},\propto \phantom{\rule{0.3em}{0ex}}\eta {\lambda}_{\left(n\right)}{E}_{0}, and \propto \phantom{\rule{0.3em}{0ex}}\eta {\lambda}_{\left(n\right)}^{2}{E}_{0}, respectively, entail the following equation:

{\left[{\phi}_{\left(n\right)}^{k}\left({\tau}_{\mathrm{ES}}\right)\right]}_{max}\phantom{\rule{2.36043pt}{0ex}}=\phantom{\rule{2.36043pt}{0ex}}a{\lambda}_{\left(n\right)}\phantom{\rule{2.36043pt}{0ex}}+\phantom{\rule{2.36043pt}{0ex}}\mathrm{a\eta \zeta}{\lambda}_{\left(n\right)}\phantom{\rule{2.36043pt}{0ex}}+\phantom{\rule{2.36043pt}{0ex}}\mathrm{a\eta}{\lambda}_{\left(n\right)}^{2}

(33)

in which *a* is an unknown scalar, and \zeta =\frac{1}{2}\sqrt{2}. For the special case *n*=0, i.e. the start of the measurements and the very cardiac cycle that the initial calibration refers to, we have:*λ*
_{(0)}=1, because by definition {\lambda}_{\left(n\right)}=\frac{{\stackrel{~}{V}}_{\left(n\right)}}{{\stackrel{~}{V}}_{\left(0\right)}} (see eq.(16)). As a result, the unknown scalar *a* cancels out if we consider the ratio \left(\right)close="">\n \n \n \n \phi \n \n \n (\n n\n )\n \n \n k\n \n \n (\n \n \n \tau \n \n \n ES\n \n \n )\n /\n \n \n \phi \n \n \n (\n 0\n )\n \n \n k\n \n \n (\n \n \n \tau \n \n \n ES\n \n \n )\n \n as function of \frac{{\stackrel{~}{V}}_{\left(n\right)}}{{\stackrel{~}{V}}_{\left(0\right)}}:

{\left[\phantom{\rule{2.36043pt}{0ex}}\frac{{\phi}_{\left(n\right)}^{k}\left({\tau}_{\mathrm{ES}}\right)}{{\phi}_{\left(0\right)}^{k}\left({\tau}_{\mathrm{ES}}\right)}\phantom{\rule{2.36043pt}{0ex}}\right]}_{max}\phantom{\rule{2.36043pt}{0ex}}=\phantom{\rule{2.36043pt}{0ex}}\frac{{\stackrel{~}{V}}_{\left(n\right)}}{{\stackrel{~}{V}}_{\left(0\right)}}\phantom{\rule{0.3em}{0ex}}\xb7\phantom{\rule{0.3em}{0ex}}\frac{1+\mathrm{\eta \zeta}+\eta \left(\frac{{\stackrel{~}{V}}_{\left(n\right)}}{{\stackrel{~}{V}}_{\left(0\right)}}\right)}{1+\mathrm{\eta \zeta}+\eta}

(34)

Using known specific conductivities of tissues from the literature, and known typical anatomical distances and sizes, the value of *η* for feature (b) has been calculated, resulting in *η*=0.035. The corresponding graph of formula (34) using *η*=0.035is the dotted green line in Figure 10. The slender area between this dotted green line and the identity line (black solid line) in Figure 10, represents the area containing all possible deviations from linearity for feature (b).

Similarly, a graph was calculated for the entire perturbation series in eq.(31), using a geometrical power series (i.e.: \sum _{\vartheta =0}^{\infty}\phantom{\rule{0.3em}{0ex}}{p}^{\vartheta}=1/(1-p) for any scalar *p*) for feature (b). This resulted in the continuous green line in Figure 10.

For each of the three features (a), (b) and (c) mentioned above, the *η* values have been calculated.

The value of *η* for the feature (a) is −0.033, and the combination of features (b) and (c), excluding feature (a) to create a worst case, yields *η*=0.142. The corresponding graphs are rendered in Figure 10 as well.

For a volume change {\stackrel{~}{V}}_{\left(n\right)} that equals two times the initial volume change {\stackrel{~}{V}}_{\left(0\right)}, the non-linear effects may cause a voltage change that is maximally 14% higher than the linear approximation. Using the linear approximation therefore leads to an *overestimation* of the stroke volume of maximally 14% in this case.

For a volume change {\stackrel{~}{V}}_{\left(n\right)} that equals half of the initial volume change {\stackrel{~}{V}}_{\left(0\right)}, the non-linear effects may cause a voltage change that is maximally 6% lower than the linear approximation. In this case, the linear approximation therefore leads to an *underestimation* of the stroke volume of maximally 6%.