The radial formulation of the bidomain model has been used to compute a reference solution, on a circular domain where the heart H occupies the region 0 < r < 1.0 and the surrounding torso T is 1.0 < r < 1.4. The reference solution was computed on a grid with a total of 2241 nodes, which gives a nodal distance h of 0.00625 mm. The time step used is At = 0.001 ms. It is clear that this resolution level would be difficult to achieve on a two- or three-dimensional domain, because the number of unknowns would become too large to handle. Table 3.3 shows convergence results for the bidomain model, using the cell model by Winslow et al. [145]. Only the errors in the transmembrane potential v are reported,and the error is measured at t = 3 ms. As above, the error is measured in the L2 norm. We see that the results are similar to those for the monodomain model, although the convergence results are slightly less favourable. For the choice 6 = 1, which combines a Godunov splitting with the backward Euler method, we see less than first-order accuracy, while the algorithm that combines Strang splitting with the Crank-Nicolson discretization gives an accuracy close to order two on the finer grids. The more variable results seen in this computation compared to the previous one may be caused by the very rapid dynamics of the applied cell model. Compared to the FitzHugh-Nagumo model, the model by Winslow et al. has a very rapid upstroke, which leads to sharp gradients in the solution. This places stricter demands on the spatial resolution required to observe consistent convergence results.

The limitations of the present approach, where a numerical solution is used as the reference for convergence estimates, were discussed above. Although the convergence of the algorithms is confirmed, our confidence in the results depends on our confidence in the correctness of the numerical reference solution. Furthermore, the tests are performed on highly idealized geometries, and we can only assume that the demonstrated convergence properties hold for more realistic experiments. Still, for a complex model like the bidomain model, it is difficult to come up with good alternatives to this type of convergence experiment. Analytical solutions can only be obtained after substantial simplifications of the equation, and the value of using these solutions for convergence tests is, therefore, also very limited. Although the results should be used with some care, the convergence tests presented here are, therefore, valuable for examining the accuracy of the complex numerical schemes that are employed.

The formulation of the equations in polar coordinates can also be performed in 3D, using the spherical coordinates (r, Q, p). By making the same assumptions as for the 2D case, i.e. spherical symmetry and scalar, constant conductivities, we obtain systems of equations very similar to those presented above. The equations can be solved as 1D problems depending only on the radius r, and used as a reference for checking the convergence of 3D solutions. However, the convergence tests in 3D are not as easy to perform as in the 2D case. To get a picture of the convergence the experiments must be performed for a number of different grid refinement levels, and in 3D this will quickly lead to a very large number of nodes, which easily exceeds the capacity of the available computers. Therefore, until computer speed and memory have increased sufficiently, it is easier to base convergence tests on 2D experiments.

3.4.2 Simulation on a 2D Slice

The results shown in previous sections were only suitable for showing the convergence properties of the numerical methods. The practical value of simulations with perfect circular symmetry is, of course, very limited. To illustrate how the numerical simulations can be used for practical purposes, we will now present a different set of experiments. The geometry used for the experiments is a two-dimensional slice through the heart and torso. The advantages of doing 2D studies, rather than 3D, are that the computational burden is much smaller and it is easier to visualize the results. The disadvantage is that the results may be unrealistic in many respects; particularly, in our case, regarding the propagation pattern. Forcing the signal propagation to occur in only one plane will not give a realistic activation pattern. However, 2D simulations can still provide insight into several important mechanisms of the heart. The goal of the present study is to investigate the relation between ischemic heart disease and shifts in the ST segment of the ECG. ST segment deflection is an important indicator for ischemic heart disease, and so the relation is important from a clinical point of view.

Geometry and data. The geometry used for the simulations is shown in Figure 3.10. The torso boundary is based upon a cross section from the Visible Human dataset [141]. This dataset consists of photographs of horizontal cross sections of a male torso body. The chosen section cuts through the middle of the ventricles. The Visible Human dataset was not used for the heart boundary, because the heart cavities are not so clear in the photographs. The boundary was traced from a different set of photographs [84].

The tissue is stimulated at three different locations, as indicated in Figure 3.10 (a). The locations represent the endings of fibres in the Purkinje network. The left ventricular sites are triggered 20 ms before the right ventricular site, the difference in timing being based on the activation map of Duirrer [33].

10 cm

Fig. 3.10. a) The heart and torso boundaries. The black dots on the endocardial surface indicate stimulation points. The plus and minus poles on the torso surface show the location where the surface potential is recorded. b) The fibre directions in the heart. The large arrows are the selected directions and the small arrows are interpolations.

The directions of the fibres are shown in Fig. 3.10b). They run tangentially to the heart surfaces. The directions are not based directly on measurements, but are selected to give a reasonable representation of the orientation in the two-dimensional slice. The out-of-plane component of the fibres is not taken into account.

3.4.3 Normal Propagation

For the simulation under consideration, there were 61167 nodes in the heart and a total of 292417 in the body. The average node density was about 0.2 mm in the heart and 0.4 mm in the body. The time step was 0.125 ms and the simulation time was 500 ms, giving a total of 4000 steps. The computational loads of the ODE and the PDE steps of the operator-splitting algorithm were approximately equal.

Figure 3.11 shows four snapshots of the simulation. After 25 ms, we see that the depolarization front has spread out from the two stimulation points on the left endocardial surface and we also see the early contribution from the third stimulation point. Notice from the position of the wave front that the propagation velocity is larger in the direction along the fibres than across them. This reflects the larger conductivity in the direction of the fibres. After 75ms, most of the tissue has been depolarized. The third image shows the situation after 220 ms. The whole heart has now depolarized and we see the beginnings of repolarization in some areas. These areas coincide with the initial depolarization. The reason that the repolarization follows the same pattern as the depolarization is that the action potential duration is the same in all locations This is not realistic, since it has been reported that the en-docardial cells, and in particular the cells near the middle of the heart wall, have longer action potential durations than epicardial cells. The effect of this is that although the depolarization of the tissue starts at the endocardium, the repolarization

Fig. 3.11. The transmembrane potential (mV) at four stages during normal propagation. (For the color version, see Figure A.6 on page 290).

wave starts at the epicardium. In our model we use identical cells throughout the myocardium, which results in both the depolarization wave and the repolarization wave starting from the endocardium. In the last snapshot the repolarization is further advanced. Notice that the repolarization front is less steep than the depolarization front. This reflects the shape of the action potential, which has a fast upstroke and a much slower repolarization.

3.4.4 Ischemia

In this section we look at the effect of changing the ionic model in parts of the tissue. Figure 3.12 shows the area where the model has been modified. The dark area represents ischemic tissue, and the ionic model is modified to reflect this. Figure 3.13 shows the normal action potential used in the healthy part of the tissue, along with the modified action potential used in the ischemic region.

Four changes have been made to the model, to reflect the changes occurring in the cells during ischemia. The extracellular potassium concentration has been increased from 4 mM to 9 mM; the maximum conductance of the fast sodium channels is reduced by 25%; the maximum conductances for Ito, IKr, and IKs have been reduced by 50%; and conductance of the time independent potassium current IKl has been reduced by 60%. These changes may not correspond exactly to the changes taking place during ischemia, but the resulting action potential is reasonable, which

Your Heart and Nutrition

Your Heart and Nutrition

Prevention is better than a cure. Learn how to cherish your heart by taking the necessary means to keep it pumping healthily and steadily through your life.

Get My Free Ebook

Post a comment