- Email: [email protected]
Programs for assessment of spatial heterogeneity of regional organ blood flow
Computer Methods and Programs in Biomedicine 55 (1998) 51 – 57
Programs for assessment of spatial heterogeneity of regional organ blood flow M. Kleen a,b,*, O. Habler a,b, B. Zwissler b, K. Messmer a a
Institute for Surgical Research, Klinikum Grosshadern, Marchioninistr. 15, 81366 Munich, Germany b Institute of Anesthesiology, Klinikum Grosshadern, Marchioninistr. 15, 81366 Munich, Germany Received 12 February 1997; received in revised form 28 July 1997; accepted 30 July 1997
Abstract Regional organ blood flow (RBF) is spatially heterogeneous. Relative dispersion (standard deviation S.D./mean) is often used to assess heterogeneity of RBF. Relative dispersion is a global measure of heterogeneity and is strongly influenced by the tissue sample size making comparisons between research groups inappropriate. Spatial correlation (SC) of blood flow is, on the other hand, averaged local self similarity. Both parameters change oppositely secondary to interventions. Fractal dimension (D) is a scale-independent measure of spatial heterogeneity and thus facilitates comparison of data. Programs for calculation of SC and D have not been published. We present two portable computer programs written in C + + for calculating SC and D. The programs were validated with six computer generated data sets of known heterogeneity. The results were in agreement with data from the literature: we conclude that the programs accurately calculate spatial correlation and fractal dimension of 1-, 2-, or 3-dimensional perfusion matrices. © 1998 Elsevier Science Ireland Ltd. Keywords: Spatial heterogeneity; Spatial correlation; Regional blood flow; Fractal dimension
1. Introduction Regional blood flow (RBF) to organs is highly heterogeneous with respect to spatial distribution [1 – 4]. Temporal heterogeneity has furthermore been demonstrated for the pulmonary microcirculation [5]. Spatial and temporal heterogeneity of * Corresponding author. Tel.: + 49 89 70954357; fax: + 49 89 70958897; e-mail: [email protected]
regional blood flow are usually quantified with the relative dispersion (RD) that is defined as standard deviation (S.D.) of blood flow divided by average blood flow. Since for calculation of RD all flow values for the organ region analyzed are averaged, RD is a global measure of heterogeneity and does not include information of local heterogeneity. The heterogeneity as assessed by RD is highly dependent on the dissection scheme that is used to divide organs into samples [3]. Higher
0169-2607/98/$19.00 © 1998 Elsevier Science Ireland Ltd. All rights reserved. PII S 0 1 6 9 - 2 6 0 7 ( 9 7 ) 0 0 0 4 7 - 3
52
M. Kleen et al. / Computer Methods and Programs in Biomedicine 55 (1998) 51–57
resolution and, thus, smaller tissue samples result in higher RD. This makes comparison of data from different research groups inappropriate. The dependence of heterogeneity on resolution, however, is the basis of the application of fractal geometry to the study of organ blood flow heterogeneity. After development by Mandelbrot [6], fractal geometry was applied to the study of organ blood flow. It was shown for many organs that heterogeneity of blood flow and the resolution of measurement of blood flow are related obeying a power law. Thus, organ perfusion may be considered a fractal process [2,7,8,3]. The fractal property of RBF may be quantified by the fractal dimension D that is defined as one minus the exponent of the power law describing the relationship of resolution and heterogeneity [2]. Since fractal dimension D, by definition, incorporates different resolutions of measurement it is independent of organ sample size and therefore facilitates comparison of data from different studies. Heterogeneity of pulmonary perfusion has, in part, been explained by the dichotomous, fractal structure of pulmonary arterial vessels [9]. On this basis, the distance between individual lung tissue samples was studied to elucidate it’s role in perfusion heterogeneity. As a result, it was shown that there is high local correlation of blood flow to neighbored pulmonary tissue samples [10]. This relation becomes weaker with distance and, at a certain distance, even negative. Spatial correlation of regional blood flow has very recently been used for assessment of heterogeneity for the myocardium [4,11]. For detailed analysis of spatial heterogeneity, spatial correlation may be calculated for all possible distances. The resulting correlation values can be plotted together with the corresponding distance and the decay of correlation with increasing distance can be analyzed [10]. Spatial correlation of directly neighbored tissue samples, however, yields a measure of averaged local perfusion heterogeneity. Conversely, relative dispersion and fractal dimension of perfusion are based on global averaging of perfusion values and are thus measures of global heterogeneity.
Changing the posture of anesthetized dogs from supine to prone decreases relative dispersion of pulmonary perfusion (i.e. decreases scale-dependent global heterogeneity) and increases fractal dimension (i.e. increases scale-independent global heterogeneity). Conversely, decreased spatial correlation of adjacent lung tissue samples indicates increased local perfusion heterogeneity [5]. We have confirmed the discrepancy between RD and spatial correlation for canine pulmonary perfusion after isovolemic hemodilution from hematocrit 36–20% (manuscript submitted to the ‘American Journal of Physiology’). From these data it can be concluded that averaged local perfusion heterogeneity (spatial correlation of adjacent tissue samples) yields additional information about perfusion heterogeneity. Therefore, for a complete analysis of regional organ perfusion and heterogeneity thereof, relative dispersion, spatial correlation, and fractal dimension should be calculated. Despite the frequently documented use for experimental studies of perfusion heterogeneity, the details of calculation of fractal dimension and spatial calculation have not been formalized in a computer program that is publicly available. Regional perfusion is generally measured with great resolution (lung: up to 1700 tissue samples per animal [5]; heart: up to 500 tissue samples per animal [1]). Therefore, the iterative processes of calculating spatial correlation and fractal dimension are extremely lengthy and virtually impossible to perform without automation. We therefore present in this paper portable C+ + programs for calculation of relative dispersion, spatial correlation, and fractal dimension of spatial (2- or 3-dimensional) or linear (1-dimensional) data.
2. Computational theory
2.1. Spatial correlation The algorithm that calculates spatial correlations for all possible distances within the bounds of the 3-dimensional arrays of blood flows implements the 3-dimensional extension of the linear
M. Kleen et al. / Computer Methods and Programs in Biomedicine 55 (1998) 51–57
correlation formula reported by Glenny in 1992 [10]: l
r(Dx, Dy, Dz)=
'
m
53
memory. Each of these data structures contains information about all flow pairs that are sepa-
n
% % % [ f(xi, yj, zk ) −F( 1] · [ f(xi + Dx, yj + Dy, zk + Dz)− F( 2] i=1 j=1 k=1 l
m
n
l
m
(1)
n
% % % [ f(xi, yj, zk ) − F( 1] · % % % [ f(xi + Dx, yj + Dy, zk + Dz)− F( 2]
i=1 j=1 k=1
2
2
i=1 j=1 k=1
The left side of Eq. (1) represents the spatial correlation for one Euclidean distance with differences in the x, y, z directions of Dx, Dy, and Dz. The actual distance is numerically defined as
Dx 2 +Dy 2 +Dz 2. The coordinate system’s dimensions are given by l, m, and n in the equation. The flow at one place within the organ is given as f(x, y, z). The paired flow located at a distance of Dx, Dy, and Dz. is f(x + Dx, y +Dy, z +Dz). F( 1 and F( 2 are the corresponding means of the flow values. i, j, and k are iteration indices. The blood flow values and coordinates are represented internally as a 3-dimensional array. This must be large enough to hold all blood flow values but, on the other hand, increasing the size of the array exponentially prolongs program execution. The central algorithm is executed 2 · n 6 times, where n is the number of array elements. Using one static array size that will always be large enough to hold all flow values is therefore unacceptable. However, the C/C + + programming language does not allow dynamic creation of 3-dimensional arrays with a size that is determined at run time. Therefore, a 1-dimensional array is used and perfusion values are stored and retrieved in a way that simulates three dimensions. The minimum size of the array is determined by analysis of the coordinates in the blood flow file. The array is then created dynamically and blood flow values are stored in it. In the first part of the algorithm, average flows F( 1 and F( 2 are calculated as follows. The blood flow array is sequentially read and for every existing flow the array is searched for one or more existing flow values at all possible distances within the bounds of the coordinate system. For every distance that occurs within the array, a separate data structure is allocated within a dynamically linked list in
rated by this distance. For each distance, averages are calculated for the flows at the beginning and at the end of the distance vector. Distances differing by more than one tenth of one unit distance are considered distinct. In a further step, the dynamic list containing all possible distances and the average flows in the organ is sequentially read and the average flows are subtracted from single flows. The same is done for the flows separated by the specific distance from the original flow. According to the enumerator of the right side of Eq. (1), the resulting differences are multiplied and summed. Simultaneously, the algorithm performs the analogous steps for calculating the square root of the summed squared differences as defined in the denominator on the right side of Eq. (1). Afterwards, the denominator and enumerator for each distance are used to calculate spatial correlations. The necessary data structures and the algorithm for calculating spatial correlations are encapsulated within one ‘class’ in a separate header file. The pseudo 3-dimensional array was also implemented as a ‘class’ and is contained in a header file. This approach facilitates re-using and maintaining the C+ + code.
2.2. Fractal dimension A process or object may display heterogeneity that is dependent on the resolution of the method of measurement that is applied. If this process is fractal this dependence is a non-linear one. The relationship of heterogeneity and resolution can be described with a simple linear equation by taking the natural logarithm (ln) of both the heterogeneity and the resolution of measurement [2]:
54
M. Kleen et al. / Computer Methods and Programs in Biomedicine 55 (1998) 51–57
Fig. 1. Spatial correlation of simulated data set from a fractal computer model for pulmonary perfusion [5]. Abscissa: distance expressed as multiples of sample size (units distance). Ordinate: spatial correlation of all simulated flow pairs that are separated by the respective distance. Modeled blood flow matrix consisted of 16 · 16 · 16 values.
ln H(r) =(1−D) · ln(r/r0) +ln H(r 0)
(2)
r and r0 are some resolution and the reference resolution. H(r) and H(r 0) are the heterogeneity measured at some resolution and at the reference resolution. D is one minus the slope of the resolution/heterogeneity plot and is called fractal dimension. The term ln H(r 0) is the intercept of this plot. It is dependent on the arbitrarily chosen reference resolution and does not affect the determination of D. Thus, fractal dimension may be determined by repeatedly measuring the heterogeneity at different resolutions, constructing a logarithmic plot of heterogeneity and resolution, and assessing the slope. The fractal dimension program accepts identical input as described for the spatial correlation program and the 3-dimensional data are represented the same way. The program calculates RD of the data at the given level of resolution. Neighbored data points in the 3-dimensional data matrix are then repetitively combined to larger units and RD is determined. Each combination to increasingly large units is performed as follows: a new 3-dimensional array
is constructed with either the x-, y-, or z-dimension divided by two or three, whichever yields an integer result. All values from the array of the preceding step are copied into the new structure while assigning the appropriate re-scaled x-, y-, or z-dimension. Two or three original values are added for each new array place, respectively. Each combination with respective RD and resolution data is stored as an element of a dynamic list. For determination of D, linear regression of ln(RD) and ln(resolution) was calculated to determine the slope of Eq. (2): for better readability, we present Eq. (2) in the standardized linear form: y=a · x+ b
(3)
where y represents ln(H(r)), a is the slope (1− D), x is ln(r/r0) and b is the intercept ln(H(r 0)). The slope is then defined as n · % xi · yi − % xi · % yi a=
n · % x − % xi 2 i
2
(4)
M. Kleen et al. / Computer Methods and Programs in Biomedicine 55 (1998) 51–57
with n being the number of data points. The intercept is calculated as % x 2i · % yi − % yi · xi · % xi
b=
n · % x 2i − % xi
(5)
2
[12]. The linear correlation coefficient is determined to assess the presence of a linear relationship and is calculated as r % xi · yi − 1/n · % xi · % yi
=
' n
n
% x 2i − 1/n · % xi
n · % yi
2
2
· % y 2i −1/
55
The necessary organization of ASCII or ANSI input files for both programs is described in the file readme.txt on the ftp-server (see Section 6). Both programs produce the specified output files with results. The format is also explained in the file readme.txt. We also provide another version of the program (fract – ap.exe) which writes only the fractal dimension and the correlation coefficient into the file d.out and appends this file with every consecutive call of the program. This facilitates collection of data when analyzing experimental series. Both programs are also provided on the ftp-server as non-portable Windows Programs which implement all described features.
4. Results from using the program (6)
[13]. From these equations, the fractal dimension D, the intercept of the plot, and the correlation coefficient are calculated. Recombinations of the original array are only included if the number of values exceeds eight in order to obtain reliable values for relative dispersion.
3. Description of programs Both programs were developed in C+ + using only language elements conforming to the ANSI standard (Watcom C/C + + compiler V. 10; Powersoft, Concord, MA). This approach enables recompilation with a C+ + compiler under any other system. Both programs were implemented as command line versions to facilitate multiple calls with DOS command batch files. The programs accept as parameters: (i) the name of the file containing blood flow data; (ii) the name of the file where output is to be written; and, in the case of the spatial correlation program, (iii) a plus ( + ) or minus (−) character to turn on or off the continuous textual output during calculations that slows down program execution. With a DOS batch file containing multiple calls to the programs with appropriate file specifications, large amounts of data can be processed without user interaction.
The fractal dimension program was validated by re-calculating relative dispersion and fractal dimension of a data set presented in full detail in the appendix of Ref. [2]. For this data set our program calculated identical fractal dimension and identical relative dispersion as reported by Glenny et al. The spatial correlation program was validated and the fractal program further tested as follows: six test data sets were produced and the programs were used to calculate spatial correlation, relative dispersion, and fractal dimension for each data set. These flow values were obtained from a fractal branching computer model of pulmonary perfusion [5]. The model incorporates gravitational forces, slight random variation of flow distribution at bifurcation sites (normal distribution of variation; mean 0.5; S.D. 0.08) and yields 3-dimensional coordinates of flows. It accurately reflects relative dispersion and spatial correlation of canine regional pulmonary perfusion. The normally distributed random number generator [14] for this model was started with a different seed for each data set. Spatial correlations of such modeled data sets were expected to display the typical pattern of spatial correlation of canine regional pulmonary perfusion: directly neighbored samples correlate closely (i.e. 0.5–0.8), and correlation gradually decreases with distance. Zero correlation is reached at about 5–15 times the
56
M. Kleen et al. / Computer Methods and Programs in Biomedicine 55 (1998) 51–57
Fig. 2. Linear plot of the natural logarithms of relative dispersion (ordinate) and relative mass of ‘samples’ (abscissa). ln(Hr ): natural logarithm of heterogeneity at some resolution. ln(r/r0): natural logarithm of resolution divided by reference resolution. Data are from a simulated data set from a fractal computer model for pulmonary perfusion [5]. The modeled blood flow matrix consisted of 16 · 16 · 16 values. The linear regression equation calculated by the fractal program is given. The slope of the plot subtracted from one is the fractal dimension. The correlation coefficient was −0.99. Note that data points for the two largest recombinations were calculated separately are not included into the linear regression analysis.
distance of adjacent tissue samples (5 – 15 units distance). The scatter of correlation values above zero (i.e. small distances) is expected to be low. For larger distances, non-homogeneous, scattered correlation values are expected. The result of an example of the validation runs of the spatial correlation program is shown in Fig. 1. As expected, the pattern of distance dependent decay of correlations were identical to the results from the original model described by Glenny as well as the data from dogs reported by Glenny et al. [10,5]. The average correlation of simulated neighboring tissue samples was 0.589 0.06 (mean9 S.D.). Glenny et al. reported spatial correlation of 0.6390.04 for datasets produced by the model [5]. Zero correlation was reached at approximately 6–15 units distance which is also not different from the values found by these authors. For greater distances, the scatter of correlations was, as expected, more than for closer distances. In Fig. 2, the results of the fractal dimension
program for the same data set are shown. The two largest recombinations shown on the right side of the plot were calculated by a modified version of the program and are not included into the linear regression analysis. The program discards recombinations with eight or less values since the calculation of relative dispersion is not reliable for very small data sets. Average fractal dimension for all data sets was 1.139 0.04 with a mean correlation coefficient of the ln(RD)/ln(r/r0) plot of− 0.989 0.01. Glenny et al. reported an average D of 1.13 (no range) of such simulated data sets [5]. From the output of both programs for simulated perfusion data, we conclude that our programs accurately calculate spatial correlation as well as fractal dimensions of 3-dimensional matrices. 1-dimensional data may also be validly assessed with the fractal dimension program since a linear time course signal data set with known fractal dimension from the literature [2] was correctly handled by our program.
M. Kleen et al. / Computer Methods and Programs in Biomedicine 55 (1998) 51–57
57
5. Hardware and software specifications
References
The compiled programs require any version of DOS, a 8086 or higher CPU. The Windows versions require MS-Windows 3.1 or higher. The source code of the non-Windows versions was written in standard (ANSI) C/C + + and can therefore be recompiled under any other operating system for that a C + + compiler is available. To the authors’ knowledge, a C+ + compiler is available for all computer systems used in biomedical research. The ‘fractal’ program occupies approximately 71 kb, the ‘spatial correlation’ program 57 kb of computer memory (compiled for IBM PC). Additionally, the data structures that are allocated during the run-time of the ‘spatial correlation’ program with a 16 · 16 · 16 data matrix require 160 kb (using a 16 bit based system). The ‘fractal’ program uses one tenth of this amount. Reading the flow values from the input file, determining array size, and initializing the flow array takes a few seconds even for very large blood flow matrices. In the ‘spatial correlation’ program, the most time consuming steps are 2 · n 6 iterations of the algorithms for assessing mean flows and differences of flows and mean flows. For a very large matrix of 16 · 16 · 16 values, program execution time is approximately 15 min on an IBM-compatible computer with a 66 MHz80486 processor. The ‘fractal’ program requires only a few seconds for data sets of this size. The Windows versions are not notably slower.
[1] B. Zwissler, R. Schosser, V. Iber, M. Weiss, C. Schwickert, P. Spengler, K. Messmer, Methodological error and spatial variability of organ blood flow measurements using radiolabeled microspheres, Res. Exp. Med. 191 (1991) 47 – 63. [2] R.W. Glenny, H.T. Robertson, S. Yamashiro, J.B. Bassingthwaighte, Applications of fractal analysis to physiology, J. Appl. Physiol. 70 (1991) 2351 – 2367. [3] J.E. McNamee, Fractal perspectives in pulmonary physiology, J. Appl. Physiol. 71 (1991) 1 – 8. [4] H. Mori, M. Chujo, S. Haruyama, H. Sakamoto, Y. Shinozaki, M. Uddin-Mohammed, A. Iida, H. Nakazawa, Local continuity of myocardial blood flow studied by monochromatic synchrotron radiation-excited X-ray fluorescence spectrometry, Circ. Res. 76 (1995) 1088 – 1100. [5] R.W. Glenny, N.L. Polissar, S. McKinney, H.T. Robertson, Temporal heterogeneity of regional pulmonary perfusion is spatially clustered, J. Appl. Physiol. 79 (1995) 986 – 1001. [6] B.B. Mandelbrot, The fractal Geometry of Nature, Freeman, New York, 1983. [7] P.O. Iversen, G. Nicolaysen, Fractals describe blood flow heterogeneity within skeletal muscle and within myocardium, Am. J. Physiol. (Heart Circ. Physiol. 37) 268 (1995) H112 – H116. [8] P.O. Iversen, G. Nicolaysen, High correlation of fractals for regional blood flows among resting and exercising skeletal muscles, Am. J. Physiol. (Heart Circ. Physiol.38) 38 (1995) H7 – H13. [9] R.W. Glenny, H.T. Robertson, Fractal properties of pulmonary blood flow: characterization of spatial heterogeneity, J. Appl. Physiol. 69 (1990) 532 – 545. [10] R.W. Glenny, Spatial correlation of regional pulmonary perfusion, J. Appl. Physiol. 72 (1992) 2378 – 2386. [11] T. Matsumoto, M. Goto, H. Tachibana, Y. Ogasawara, K. Tsujioka, F. Kajiya, Microheterogeneity of myocardial blood flow in rabbit hearts during normoxic and hypoxic states, Am. J. Physiol. (Heart Circ. Physiol. 39) 270 (1996) H435 – H441. [12] B. Ramm, G. Hofmann (Eds.), Korrelation, in: Biomathematik, Ferdinand Enke Verlag, Stuttgart, 1982, pp. 75 – 77. [13] B. Ramm, G. Hofmann (Eds.), Regressionsgerade, in: Biomathematik, Ferdinand Enke Verlag, Stuttgart, 1982, pp. 78 – 79. [14] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling (Eds.), Random numbers, in: Numerical Recipes in Pascal, Cambridge University Press, Cambridge, 1992, pp. 212 – 253.
6. Mode of availability By the time of publication, the executable files and complete source codes are available on the anonymous FTP-server of the Leibniz – Rechenzentrum in Munich, Germany: ftp:// ftp.lrz/muenchen.de/pub/science/medicine/surgicalresearch/heterogeneity.
.
Programs for assessment of spatial heterogeneity of regional organ blood flow M. Kleen a,b,*, O. Habler a,b, B. Zwissler b, K. Messmer a a
Institute for Surgical Research, Klinikum Grosshadern, Marchioninistr. 15, 81366 Munich, Germany b Institute of Anesthesiology, Klinikum Grosshadern, Marchioninistr. 15, 81366 Munich, Germany Received 12 February 1997; received in revised form 28 July 1997; accepted 30 July 1997
Abstract Regional organ blood flow (RBF) is spatially heterogeneous. Relative dispersion (standard deviation S.D./mean) is often used to assess heterogeneity of RBF. Relative dispersion is a global measure of heterogeneity and is strongly influenced by the tissue sample size making comparisons between research groups inappropriate. Spatial correlation (SC) of blood flow is, on the other hand, averaged local self similarity. Both parameters change oppositely secondary to interventions. Fractal dimension (D) is a scale-independent measure of spatial heterogeneity and thus facilitates comparison of data. Programs for calculation of SC and D have not been published. We present two portable computer programs written in C + + for calculating SC and D. The programs were validated with six computer generated data sets of known heterogeneity. The results were in agreement with data from the literature: we conclude that the programs accurately calculate spatial correlation and fractal dimension of 1-, 2-, or 3-dimensional perfusion matrices. © 1998 Elsevier Science Ireland Ltd. Keywords: Spatial heterogeneity; Spatial correlation; Regional blood flow; Fractal dimension
1. Introduction Regional blood flow (RBF) to organs is highly heterogeneous with respect to spatial distribution [1 – 4]. Temporal heterogeneity has furthermore been demonstrated for the pulmonary microcirculation [5]. Spatial and temporal heterogeneity of * Corresponding author. Tel.: + 49 89 70954357; fax: + 49 89 70958897; e-mail: [email protected]
regional blood flow are usually quantified with the relative dispersion (RD) that is defined as standard deviation (S.D.) of blood flow divided by average blood flow. Since for calculation of RD all flow values for the organ region analyzed are averaged, RD is a global measure of heterogeneity and does not include information of local heterogeneity. The heterogeneity as assessed by RD is highly dependent on the dissection scheme that is used to divide organs into samples [3]. Higher
0169-2607/98/$19.00 © 1998 Elsevier Science Ireland Ltd. All rights reserved. PII S 0 1 6 9 - 2 6 0 7 ( 9 7 ) 0 0 0 4 7 - 3
52
M. Kleen et al. / Computer Methods and Programs in Biomedicine 55 (1998) 51–57
resolution and, thus, smaller tissue samples result in higher RD. This makes comparison of data from different research groups inappropriate. The dependence of heterogeneity on resolution, however, is the basis of the application of fractal geometry to the study of organ blood flow heterogeneity. After development by Mandelbrot [6], fractal geometry was applied to the study of organ blood flow. It was shown for many organs that heterogeneity of blood flow and the resolution of measurement of blood flow are related obeying a power law. Thus, organ perfusion may be considered a fractal process [2,7,8,3]. The fractal property of RBF may be quantified by the fractal dimension D that is defined as one minus the exponent of the power law describing the relationship of resolution and heterogeneity [2]. Since fractal dimension D, by definition, incorporates different resolutions of measurement it is independent of organ sample size and therefore facilitates comparison of data from different studies. Heterogeneity of pulmonary perfusion has, in part, been explained by the dichotomous, fractal structure of pulmonary arterial vessels [9]. On this basis, the distance between individual lung tissue samples was studied to elucidate it’s role in perfusion heterogeneity. As a result, it was shown that there is high local correlation of blood flow to neighbored pulmonary tissue samples [10]. This relation becomes weaker with distance and, at a certain distance, even negative. Spatial correlation of regional blood flow has very recently been used for assessment of heterogeneity for the myocardium [4,11]. For detailed analysis of spatial heterogeneity, spatial correlation may be calculated for all possible distances. The resulting correlation values can be plotted together with the corresponding distance and the decay of correlation with increasing distance can be analyzed [10]. Spatial correlation of directly neighbored tissue samples, however, yields a measure of averaged local perfusion heterogeneity. Conversely, relative dispersion and fractal dimension of perfusion are based on global averaging of perfusion values and are thus measures of global heterogeneity.
Changing the posture of anesthetized dogs from supine to prone decreases relative dispersion of pulmonary perfusion (i.e. decreases scale-dependent global heterogeneity) and increases fractal dimension (i.e. increases scale-independent global heterogeneity). Conversely, decreased spatial correlation of adjacent lung tissue samples indicates increased local perfusion heterogeneity [5]. We have confirmed the discrepancy between RD and spatial correlation for canine pulmonary perfusion after isovolemic hemodilution from hematocrit 36–20% (manuscript submitted to the ‘American Journal of Physiology’). From these data it can be concluded that averaged local perfusion heterogeneity (spatial correlation of adjacent tissue samples) yields additional information about perfusion heterogeneity. Therefore, for a complete analysis of regional organ perfusion and heterogeneity thereof, relative dispersion, spatial correlation, and fractal dimension should be calculated. Despite the frequently documented use for experimental studies of perfusion heterogeneity, the details of calculation of fractal dimension and spatial calculation have not been formalized in a computer program that is publicly available. Regional perfusion is generally measured with great resolution (lung: up to 1700 tissue samples per animal [5]; heart: up to 500 tissue samples per animal [1]). Therefore, the iterative processes of calculating spatial correlation and fractal dimension are extremely lengthy and virtually impossible to perform without automation. We therefore present in this paper portable C+ + programs for calculation of relative dispersion, spatial correlation, and fractal dimension of spatial (2- or 3-dimensional) or linear (1-dimensional) data.
2. Computational theory
2.1. Spatial correlation The algorithm that calculates spatial correlations for all possible distances within the bounds of the 3-dimensional arrays of blood flows implements the 3-dimensional extension of the linear
M. Kleen et al. / Computer Methods and Programs in Biomedicine 55 (1998) 51–57
correlation formula reported by Glenny in 1992 [10]: l
r(Dx, Dy, Dz)=
'
m
53
memory. Each of these data structures contains information about all flow pairs that are sepa-
n
% % % [ f(xi, yj, zk ) −F( 1] · [ f(xi + Dx, yj + Dy, zk + Dz)− F( 2] i=1 j=1 k=1 l
m
n
l
m
(1)
n
% % % [ f(xi, yj, zk ) − F( 1] · % % % [ f(xi + Dx, yj + Dy, zk + Dz)− F( 2]
i=1 j=1 k=1
2
2
i=1 j=1 k=1
The left side of Eq. (1) represents the spatial correlation for one Euclidean distance with differences in the x, y, z directions of Dx, Dy, and Dz. The actual distance is numerically defined as
Dx 2 +Dy 2 +Dz 2. The coordinate system’s dimensions are given by l, m, and n in the equation. The flow at one place within the organ is given as f(x, y, z). The paired flow located at a distance of Dx, Dy, and Dz. is f(x + Dx, y +Dy, z +Dz). F( 1 and F( 2 are the corresponding means of the flow values. i, j, and k are iteration indices. The blood flow values and coordinates are represented internally as a 3-dimensional array. This must be large enough to hold all blood flow values but, on the other hand, increasing the size of the array exponentially prolongs program execution. The central algorithm is executed 2 · n 6 times, where n is the number of array elements. Using one static array size that will always be large enough to hold all flow values is therefore unacceptable. However, the C/C + + programming language does not allow dynamic creation of 3-dimensional arrays with a size that is determined at run time. Therefore, a 1-dimensional array is used and perfusion values are stored and retrieved in a way that simulates three dimensions. The minimum size of the array is determined by analysis of the coordinates in the blood flow file. The array is then created dynamically and blood flow values are stored in it. In the first part of the algorithm, average flows F( 1 and F( 2 are calculated as follows. The blood flow array is sequentially read and for every existing flow the array is searched for one or more existing flow values at all possible distances within the bounds of the coordinate system. For every distance that occurs within the array, a separate data structure is allocated within a dynamically linked list in
rated by this distance. For each distance, averages are calculated for the flows at the beginning and at the end of the distance vector. Distances differing by more than one tenth of one unit distance are considered distinct. In a further step, the dynamic list containing all possible distances and the average flows in the organ is sequentially read and the average flows are subtracted from single flows. The same is done for the flows separated by the specific distance from the original flow. According to the enumerator of the right side of Eq. (1), the resulting differences are multiplied and summed. Simultaneously, the algorithm performs the analogous steps for calculating the square root of the summed squared differences as defined in the denominator on the right side of Eq. (1). Afterwards, the denominator and enumerator for each distance are used to calculate spatial correlations. The necessary data structures and the algorithm for calculating spatial correlations are encapsulated within one ‘class’ in a separate header file. The pseudo 3-dimensional array was also implemented as a ‘class’ and is contained in a header file. This approach facilitates re-using and maintaining the C+ + code.
2.2. Fractal dimension A process or object may display heterogeneity that is dependent on the resolution of the method of measurement that is applied. If this process is fractal this dependence is a non-linear one. The relationship of heterogeneity and resolution can be described with a simple linear equation by taking the natural logarithm (ln) of both the heterogeneity and the resolution of measurement [2]:
54
M. Kleen et al. / Computer Methods and Programs in Biomedicine 55 (1998) 51–57
Fig. 1. Spatial correlation of simulated data set from a fractal computer model for pulmonary perfusion [5]. Abscissa: distance expressed as multiples of sample size (units distance). Ordinate: spatial correlation of all simulated flow pairs that are separated by the respective distance. Modeled blood flow matrix consisted of 16 · 16 · 16 values.
ln H(r) =(1−D) · ln(r/r0) +ln H(r 0)
(2)
r and r0 are some resolution and the reference resolution. H(r) and H(r 0) are the heterogeneity measured at some resolution and at the reference resolution. D is one minus the slope of the resolution/heterogeneity plot and is called fractal dimension. The term ln H(r 0) is the intercept of this plot. It is dependent on the arbitrarily chosen reference resolution and does not affect the determination of D. Thus, fractal dimension may be determined by repeatedly measuring the heterogeneity at different resolutions, constructing a logarithmic plot of heterogeneity and resolution, and assessing the slope. The fractal dimension program accepts identical input as described for the spatial correlation program and the 3-dimensional data are represented the same way. The program calculates RD of the data at the given level of resolution. Neighbored data points in the 3-dimensional data matrix are then repetitively combined to larger units and RD is determined. Each combination to increasingly large units is performed as follows: a new 3-dimensional array
is constructed with either the x-, y-, or z-dimension divided by two or three, whichever yields an integer result. All values from the array of the preceding step are copied into the new structure while assigning the appropriate re-scaled x-, y-, or z-dimension. Two or three original values are added for each new array place, respectively. Each combination with respective RD and resolution data is stored as an element of a dynamic list. For determination of D, linear regression of ln(RD) and ln(resolution) was calculated to determine the slope of Eq. (2): for better readability, we present Eq. (2) in the standardized linear form: y=a · x+ b
(3)
where y represents ln(H(r)), a is the slope (1− D), x is ln(r/r0) and b is the intercept ln(H(r 0)). The slope is then defined as n · % xi · yi − % xi · % yi a=
n · % x − % xi 2 i
2
(4)
M. Kleen et al. / Computer Methods and Programs in Biomedicine 55 (1998) 51–57
with n being the number of data points. The intercept is calculated as % x 2i · % yi − % yi · xi · % xi
b=
n · % x 2i − % xi
(5)
2
[12]. The linear correlation coefficient is determined to assess the presence of a linear relationship and is calculated as r % xi · yi − 1/n · % xi · % yi
=
' n
n
% x 2i − 1/n · % xi
n · % yi
2
2
· % y 2i −1/
55
The necessary organization of ASCII or ANSI input files for both programs is described in the file readme.txt on the ftp-server (see Section 6). Both programs produce the specified output files with results. The format is also explained in the file readme.txt. We also provide another version of the program (fract – ap.exe) which writes only the fractal dimension and the correlation coefficient into the file d.out and appends this file with every consecutive call of the program. This facilitates collection of data when analyzing experimental series. Both programs are also provided on the ftp-server as non-portable Windows Programs which implement all described features.
4. Results from using the program (6)
[13]. From these equations, the fractal dimension D, the intercept of the plot, and the correlation coefficient are calculated. Recombinations of the original array are only included if the number of values exceeds eight in order to obtain reliable values for relative dispersion.
3. Description of programs Both programs were developed in C+ + using only language elements conforming to the ANSI standard (Watcom C/C + + compiler V. 10; Powersoft, Concord, MA). This approach enables recompilation with a C+ + compiler under any other system. Both programs were implemented as command line versions to facilitate multiple calls with DOS command batch files. The programs accept as parameters: (i) the name of the file containing blood flow data; (ii) the name of the file where output is to be written; and, in the case of the spatial correlation program, (iii) a plus ( + ) or minus (−) character to turn on or off the continuous textual output during calculations that slows down program execution. With a DOS batch file containing multiple calls to the programs with appropriate file specifications, large amounts of data can be processed without user interaction.
The fractal dimension program was validated by re-calculating relative dispersion and fractal dimension of a data set presented in full detail in the appendix of Ref. [2]. For this data set our program calculated identical fractal dimension and identical relative dispersion as reported by Glenny et al. The spatial correlation program was validated and the fractal program further tested as follows: six test data sets were produced and the programs were used to calculate spatial correlation, relative dispersion, and fractal dimension for each data set. These flow values were obtained from a fractal branching computer model of pulmonary perfusion [5]. The model incorporates gravitational forces, slight random variation of flow distribution at bifurcation sites (normal distribution of variation; mean 0.5; S.D. 0.08) and yields 3-dimensional coordinates of flows. It accurately reflects relative dispersion and spatial correlation of canine regional pulmonary perfusion. The normally distributed random number generator [14] for this model was started with a different seed for each data set. Spatial correlations of such modeled data sets were expected to display the typical pattern of spatial correlation of canine regional pulmonary perfusion: directly neighbored samples correlate closely (i.e. 0.5–0.8), and correlation gradually decreases with distance. Zero correlation is reached at about 5–15 times the
56
M. Kleen et al. / Computer Methods and Programs in Biomedicine 55 (1998) 51–57
Fig. 2. Linear plot of the natural logarithms of relative dispersion (ordinate) and relative mass of ‘samples’ (abscissa). ln(Hr ): natural logarithm of heterogeneity at some resolution. ln(r/r0): natural logarithm of resolution divided by reference resolution. Data are from a simulated data set from a fractal computer model for pulmonary perfusion [5]. The modeled blood flow matrix consisted of 16 · 16 · 16 values. The linear regression equation calculated by the fractal program is given. The slope of the plot subtracted from one is the fractal dimension. The correlation coefficient was −0.99. Note that data points for the two largest recombinations were calculated separately are not included into the linear regression analysis.
distance of adjacent tissue samples (5 – 15 units distance). The scatter of correlation values above zero (i.e. small distances) is expected to be low. For larger distances, non-homogeneous, scattered correlation values are expected. The result of an example of the validation runs of the spatial correlation program is shown in Fig. 1. As expected, the pattern of distance dependent decay of correlations were identical to the results from the original model described by Glenny as well as the data from dogs reported by Glenny et al. [10,5]. The average correlation of simulated neighboring tissue samples was 0.589 0.06 (mean9 S.D.). Glenny et al. reported spatial correlation of 0.6390.04 for datasets produced by the model [5]. Zero correlation was reached at approximately 6–15 units distance which is also not different from the values found by these authors. For greater distances, the scatter of correlations was, as expected, more than for closer distances. In Fig. 2, the results of the fractal dimension
program for the same data set are shown. The two largest recombinations shown on the right side of the plot were calculated by a modified version of the program and are not included into the linear regression analysis. The program discards recombinations with eight or less values since the calculation of relative dispersion is not reliable for very small data sets. Average fractal dimension for all data sets was 1.139 0.04 with a mean correlation coefficient of the ln(RD)/ln(r/r0) plot of− 0.989 0.01. Glenny et al. reported an average D of 1.13 (no range) of such simulated data sets [5]. From the output of both programs for simulated perfusion data, we conclude that our programs accurately calculate spatial correlation as well as fractal dimensions of 3-dimensional matrices. 1-dimensional data may also be validly assessed with the fractal dimension program since a linear time course signal data set with known fractal dimension from the literature [2] was correctly handled by our program.
M. Kleen et al. / Computer Methods and Programs in Biomedicine 55 (1998) 51–57
57
5. Hardware and software specifications
References
The compiled programs require any version of DOS, a 8086 or higher CPU. The Windows versions require MS-Windows 3.1 or higher. The source code of the non-Windows versions was written in standard (ANSI) C/C + + and can therefore be recompiled under any other operating system for that a C + + compiler is available. To the authors’ knowledge, a C+ + compiler is available for all computer systems used in biomedical research. The ‘fractal’ program occupies approximately 71 kb, the ‘spatial correlation’ program 57 kb of computer memory (compiled for IBM PC). Additionally, the data structures that are allocated during the run-time of the ‘spatial correlation’ program with a 16 · 16 · 16 data matrix require 160 kb (using a 16 bit based system). The ‘fractal’ program uses one tenth of this amount. Reading the flow values from the input file, determining array size, and initializing the flow array takes a few seconds even for very large blood flow matrices. In the ‘spatial correlation’ program, the most time consuming steps are 2 · n 6 iterations of the algorithms for assessing mean flows and differences of flows and mean flows. For a very large matrix of 16 · 16 · 16 values, program execution time is approximately 15 min on an IBM-compatible computer with a 66 MHz80486 processor. The ‘fractal’ program requires only a few seconds for data sets of this size. The Windows versions are not notably slower.
[1] B. Zwissler, R. Schosser, V. Iber, M. Weiss, C. Schwickert, P. Spengler, K. Messmer, Methodological error and spatial variability of organ blood flow measurements using radiolabeled microspheres, Res. Exp. Med. 191 (1991) 47 – 63. [2] R.W. Glenny, H.T. Robertson, S. Yamashiro, J.B. Bassingthwaighte, Applications of fractal analysis to physiology, J. Appl. Physiol. 70 (1991) 2351 – 2367. [3] J.E. McNamee, Fractal perspectives in pulmonary physiology, J. Appl. Physiol. 71 (1991) 1 – 8. [4] H. Mori, M. Chujo, S. Haruyama, H. Sakamoto, Y. Shinozaki, M. Uddin-Mohammed, A. Iida, H. Nakazawa, Local continuity of myocardial blood flow studied by monochromatic synchrotron radiation-excited X-ray fluorescence spectrometry, Circ. Res. 76 (1995) 1088 – 1100. [5] R.W. Glenny, N.L. Polissar, S. McKinney, H.T. Robertson, Temporal heterogeneity of regional pulmonary perfusion is spatially clustered, J. Appl. Physiol. 79 (1995) 986 – 1001. [6] B.B. Mandelbrot, The fractal Geometry of Nature, Freeman, New York, 1983. [7] P.O. Iversen, G. Nicolaysen, Fractals describe blood flow heterogeneity within skeletal muscle and within myocardium, Am. J. Physiol. (Heart Circ. Physiol. 37) 268 (1995) H112 – H116. [8] P.O. Iversen, G. Nicolaysen, High correlation of fractals for regional blood flows among resting and exercising skeletal muscles, Am. J. Physiol. (Heart Circ. Physiol.38) 38 (1995) H7 – H13. [9] R.W. Glenny, H.T. Robertson, Fractal properties of pulmonary blood flow: characterization of spatial heterogeneity, J. Appl. Physiol. 69 (1990) 532 – 545. [10] R.W. Glenny, Spatial correlation of regional pulmonary perfusion, J. Appl. Physiol. 72 (1992) 2378 – 2386. [11] T. Matsumoto, M. Goto, H. Tachibana, Y. Ogasawara, K. Tsujioka, F. Kajiya, Microheterogeneity of myocardial blood flow in rabbit hearts during normoxic and hypoxic states, Am. J. Physiol. (Heart Circ. Physiol. 39) 270 (1996) H435 – H441. [12] B. Ramm, G. Hofmann (Eds.), Korrelation, in: Biomathematik, Ferdinand Enke Verlag, Stuttgart, 1982, pp. 75 – 77. [13] B. Ramm, G. Hofmann (Eds.), Regressionsgerade, in: Biomathematik, Ferdinand Enke Verlag, Stuttgart, 1982, pp. 78 – 79. [14] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling (Eds.), Random numbers, in: Numerical Recipes in Pascal, Cambridge University Press, Cambridge, 1992, pp. 212 – 253.
6. Mode of availability By the time of publication, the executable files and complete source codes are available on the anonymous FTP-server of the Leibniz – Rechenzentrum in Munich, Germany: ftp:// ftp.lrz/muenchen.de/pub/science/medicine/surgicalresearch/heterogeneity.
.