1 Introduction

The ultrastructure of cardiac myocytes changes under pathological conditions such as dilated cardiomyopathy, diabetes and myocardial infarctions. The fission-fusion process of mitochondria is disrupted and they can become smaller and more numerous (???), although giant mitochondria can also form (?). Myofibrils can become disarranged and degenerated in diabetes (?), and their function is impaired due to the proliferation of mitochondria . The nuclei can also grow and change shape (?). Myocardial infarction causes the TATS to become disorganised and disrupted (??), while in dilated cardiomyopathy it has been shown to proliferate (?).

The relationship between cardiac cell structure and function is not yet fully understood. One way of exploring this relationship is the use of anatomically accurate computational models, which can be used to simulate and analyse cardiac function for different geometeries. Although extensive work has been done on anatomically accurate computational models of the heart at tissue and organ level (?), there are no such models of cardiac cells. We aim to create a three-dimensional computational model of a cardiac myocyte, incorporating realistic geometry.

Cardiac myocytes are all unique, but they have common features. These features have been qualitatively described in the literature, and scalar factors such as volume fractions and individual organelle dimensions have been quantified (see Section 3.1, but little has yet been done to quantify the arrangement of the organelles within the cell and noone has yet attempted to simulate it. We will use stochastic techniques to quantify and simulate this arrangment.

In the present study we have focused on the arrangement of myofibrils and mitochondria in crosssection. We used spatial point patterns to analyse and simulate the arrangement of organelle centroids.

2 Methods

2.1 Tissue processing and data acquisition

Adult male Wistar rats (200-250 g body weight) were used with ethics approval. The left ventricle was injected with 0.2 mL heparin / 0.4 mL saline to stop the heart, which was then extracted and placed in ice-cold saline. The ascending aorta was isolated, cannulated and connected to a Langendorff apparatus operated at 60 cm hydraulic pressure. The heart was perfused with Krebb’s solution to re-initiate contraction, then with St Thomas’ solution to stop the heart in systole. Next it was perfused with fixative (2.5% gluteraldehyde, 2% formaldehyde, 50 mM CaCl_{2} , in 0.15 M sodium cacodylate buffer), then the ventricle walls were isolated, cut into small pieces incubated in vials of fixative. The tissue was incubate in 50 mM CaCl_{2} in 0.15 M sodium cacodylate buffer overnight on ice. The next day, the tissue was incubated in osmium ferrocyanide for 4 hrs, washed in 0.15 M sodium cacodylate buffer, stained with 2% UA (aqueous) rinsed in DDW then dehydrate in a series of alcohols before being brought up to room temperature in acetone. The tissue was then infiltrated with epoxy resin and embedded.

Thin sections (70-80 nm) were cut with and poststained with 2% UA (aqueous) and lead citrate. The sections were examined under a Tecnai Spirit Transmission Electron Microscope from FEI and images were taken using Olympus-Soft Imaging Systems Morada. Images were taken at a magnification of × 13,500 and stitched together in Gimp to form a single image of a cell in cross-section. Organelle and cell boundaries were manually digitised in Zinc and used to calculate area, diameter (largest distance between two boundary points) and centroids. The centroid patterns were imported into R for analysis as spatial point patterns.

2.2 Spatial point pattern statistics and models

A spatial point process is a random pattern of points in d-dimensional space [ref?]. These points can have information associated with them – called a mark. Marks can either be numeric (for example the radius or area of object the point represents) or categorical (to differentiate between objects of different types). Patterns with categorical marks are called multi-type point patterns (Figure 1). In this study, the organelle arrangement was treated as a multi-type point pattern, where points were marked either ’myofibril’ or ’mitochondrion’.

The arrangement of points is characterised by the trends in the intensity of points (how many points there are within a given area) and by the interaction between points (how the positioning of one point affects the positioning of other points around it) (??).

This section gives a brief description of the functional aspects of spatial point pattern statistics and modelling techniques. ?, ? and ? provided the basis for this section; refer to them for more detailed explanations of the theory, for further examples of spatial point pattern analyses and for information on how to conduct such analyses in R.

Figure 1: An example of a multi-type point pattern, with categorical marks of two types.

In this study, each cell cross-section was treated as a point pattern. The sarcolemma defined the window of observation and the organelle centroids were treated as localised points, with categorical marks attached to differentiate between myofibrils and mitochondria. The patterns were analysed and modelled in R using the spatstat library (?).

2.2.1 Characterising Intensity

The intensity (local point density) of points can vary between different areas in a point pattern. Even patterns generated according to a homogeneous distribution (as the pattern in Figure ?? was) show small amounts of variation due to the stochastic nature of point patterns. The intensity can also vary according to an underlying trend (called the intensity function) – for example the pattern in Figure ??, in which the intensity function is a quadratic function of the x coordinate.

(a)(b)

(c)(d)

Figure 2: Example patterns: (a) a pattern generated according to a homogeneous intensity distribution and

(b) a pattern with an intensity function that is quadratic in x). Scale bars indicate the intensity in

points per unit^{2}; windows are both size 1 unit x 1 unit.

(c) and (d) represent the intensity functions in x for patterns (a) and (b) respectively.

The quadrat test is a common method for determining whether the distribution of points can be considered homogeneous (constant) throughout the window of observation. However in this study it did not give any useful results, due to the relatively small number of points in each pattern [ref?].

Intensity can also be characterised by fitting intensity models to the underlying distribution. The Akaike Information Criterion (AIC) can be used to determine the simplest model that best fits the data. Parsimonious models with a low error give low AIC values. These values are not absolute, so they can only be used to compare models fitted to the same pattern. To illustrate, a polynomial (quadratic) intensity models fitted to the pattern in Figure ?? is shown in Figure 4. The results for models fitted to the cell data are given in Section 3.2.1.

2.2.2 Characteristing Interaction

The interaction between points can be characterised as attraction (points cluster together) and repulsion (points form a regular pattern). One type of regularity is the hardcore process – when the points must lie at least r_{min} apart. If there is no interaction, points do not affect each other, so fall according to a Poisson distribution. The different types of interaction are illustrated in Figure 3.

The interaction can be investigated using the cumulative distribution function for the nearest neighbour distances ( 1) . This is called the Gfunction and denoted G(r). The nearest neighbour distances are the Euclidean distances between each point and its nearest neighbour.

(a) (b) (c)

Figure 3: Point patterns with an overall homogeneous intensity, illustrating the different types of interaction between points: (a) attraction (clustering), (b) repulsion (regularity) and (c) no interaction (‘Poisson’).

Figure 4: A polynomial (quadratic) intensity models fitted to the pattern in Figure ??. The AIC value

(AIC = 29.5) is much less than that of the homogeneous model (AIC = 140.3), indicating that it

is a more appropriate model.

G(r) = \frac{1}{N} ∑_{i=1}^{N} 1 (d(x_{i}) ≤ r) for r ≥ 0 (1)

where d(x_{i}) is the distance of point x_{i} to its nearest neighbour, and 1(expression) is the indicator function, which has the value 1 if the expression is true and 0 if false. When G(r) is above the expected Poisson curve for at small r values, it indicates clustering (more points closer together) and when it lies below the expected curve it indicates regularity (fewer points close together). Figure 5 shows the G-functions for the example patterns shown in Figures ??, ?? and ??. Figure ?? shows that the pattern in Figure ?? is likely to come from a hardcore process, as there are no points within r = 0.18 of each other.

Pairwise distance functions can be used to characterise the interaction for points beyond the nearest neighbours. Ripley’s K -function , denoted K(r), is the cumulative distance function of the distances between all pairs of points, normalised by the local intensity function.

K(r) = \frac{1}{λN} ∑_{i}^{N} ∑_{j,j≠i}^{N} 1(d_{ij} ≤ r)  (2)

where λ̂ is the estimated intensity, d_{ij} is the distance between points i and j. For a homogeneous Poisson distribution, the theoretical value for the K -function is πr^{2}.

The pair correlation function (PCF), denoted g(r), is the derivative of the K -function normalised by 2πr (expected curve from a Poisson process).

g(r) = \frac{1}{2πr} \frac{dK}{dr}  (3)

The PCF indicates how many pairs of points lie at a distance r apart, proportional to the expected number for a Poisson process. Thus when r is small, PCF values well above 1 indicate clustering (points are more likely to lie close together), and values well below one indicate regularity (points are less likely to lie close together). Figure 5 shows the PCFs for the patterns in Figures ??, ?? and ??.

For multi-type marked point patterns, the G-function and PCF can be calculated seperately for distances between points of the same type, and also as cross-functions – for distances between points of different types. The cross functions are denoted as G_{ij}(r), K_{ij}(r), and g_{ij}(r). G_{ij}(r) is the cumulative distribution function of the distance between each point of one type (i) and the point of the other type (j ) closest to it (or vice versa). K_{ij}(r) is the average number of type j points within r of a point of type i, normalised by the local intensity of type j points(λ_{j}) – or vice versa. The cross-PCF, g_{ij}(r) is calculated from K_{ij}(r), as described in (3).

In this study, the G-function and PCF were calculated separately for the distances between myofibrils and the distances between mitochondria. The corresponding cross-functions were also calculated, for myofibril-mitochondrion distances. The results are shown in Figure 9 in Section 3.2.1.

Due to the stochastic nature of point processes, it is helpful to compare observed functions to a range of acceptable values rather than a single theoretical curve. This can be done using Monte Carlo envelopes – found by generating 99 simulations of a homogeneous Poisson process, calculating the corresponding G-function or PCF curves, and taking the highest and lowest function values for each value of r. If the function for the observed pattern does not lie within the envelopes, the model is not considered reasonable, as it does not capture all the information in the data. Monte Carlo tests can also be used to determine whether or not

(a)(b)(c)

(d)(e)(f)

Figure 5: G-function and PCF for the example patterns previously shown in Figures ??, ?? and ??, showing respectively (a) & (d) clustering (obs above theo), (b) & (e) no interaction (obs close to theo) and (c) & (f) regularity (obs = 0 until r ≈ 0.18, so no points lie within that distance of each other). The curve ‘obs’ is the observed function for the pattern, and the dotted line ‘theo’ is the theoretical Poisson curve (for the G-function,

G_{theo}(r) = 1 − e^{−λπr^{2}} and for the PCF, g_{theo}(r) = 1).

an observed pattern is consistent with a particular model, by calculating Monte Carlo envelopes from 99 simulations of that model.

2.2.3 Modelling interaction

Organelles have a finite radius, so the centroids cannot come within a certain distance r_{min} of each other. This means they must be modelled with a hardcore process. A simple hard-core model is the thinned Poisson process. Points are generated according to the Poisson process (no interaction), but if they fall within the hardcore radius (r_{min}) of another point, they are removed.

In this study, a thinned Poisson model was fitted to the myofibril and mitochondrion centroid patterns from each cell. The average intensity and minimum distance between points were set to be identical to the observed centroid pattern. For each pattern, the model was compared to the observed pattern using a Monte Carlo test with the G-function and PCF.

The thinned Poisson process does not involve any interaction beyond the hardcore radius. To model a region of futher interaction, a Strauss-hardcore Gibbs model can be used. Gibbs models are governed by the intensity function and the pair potential function, ϕ(r), which defines the interaction. For a Strauss-Hardcore process:

ϕ(r) = \{ \begin{matrix}∞ & r ≤ r_{hc} \\ β & r_{hc} < r ≤ r_{S} \\ 0 & r > r_{S} \end{matrix}

A plot of ϕ(r) for the Strauss-hardcore process is shown in Figure 6. This pair potential function means that below a certain radius there is infinite repulsion, so a point cannot fall within a hardcore radius r_{hc} of another. The range r_{hc} < r ≤ r_{S}

Figure 6: Pair potential function for a Strauss-hardcore model.

defines the region of Strauss interaction. Points within this distance of each other will experience repulsion, the strength of which is defined by β. In theory, β can vary within this range, however in this study it was assumed to be constant. For a more detailed explanation of the theory behind Gibbs models see Chapter 20 in ? and Section 6.6 in ?.

When fitting a Strauss-hardcore model in spatstat, the hardcore radius and the Strauss interaction radius are user-defined, and only the intensity distribution and Strauss repulsion strength (β) are fitted to the data. Appropriate values for the hardrcore and Strauss radii were found by varying both and analysing the resulting model error for each combination of values (using the Monte Carlo envelopes). It was soon established that a hardcore radius of r_{hc} = 0.98× r_{min} (r_{min} = minimum distance between two points in the pattern) was best able to reproduce the observed hardcore behaviour, so only the Strauss interaction radius had to be optimised for each pattern.

Strauss-hardcore models were first fitted to the myofibril and mitochondrion separately and then together as a multi-Strauss-hardcore process. The latter required finding separate Strauss radius values for myofibril-myofibril, mitochondrionmitochondrion and myofibril-mitochondrion interactions. Hardcore radii were set to 0.98 × r_{min}.

Table 1: Average cell and organelle size statistics across all the cells: cell area, density of centroids, myofibril/mitochondrion percentage of the cell area (analogous to volume fraction), and average individual organelle diameters and areas. Standard deviations are given as a percentage of the mean. Corresponding

statistics found in the literature are shown in grey.

Cell area (µm)  Organelle Density (/µm^{2})  % of Cell area  Diameter ( µm)  Area ( µm^{2})
Med  Mean %SD±Med Mean %SD±
63.5±45% myofibrils 1.42±19% 50.4 ±7%  1.02  1.11 43%±0.30 0.37 69%±
-  54-61^{a}  0.7-1.1^{b}  ±  -
mitochondria 1.92 ±20% 33.9±10% 0.68  0.73 48%±0.15 0.19 81%±
approx. 1^{c} 26-36^{d}  0.5-1^{e}  ± 0.37,0.57^{f}

a???????

b???

c?

d????????

e?

f? and ? respectively

3 Results and Discussion

3.1 Size data

The average statistics for cell and organelle sizes are shown in Table 1, together with analogous results described in the literature. In general the values we obtained were comparable to previously recorded values. Values in the literature vary considerably between studies, and this variatiation probably accounts for the minor discrepencies. We took this to mean that we were correctly identifying organelle boundaries, so our centroid data could be considered reasonably accurate. When comparing individual organelle sizes from all the cells, there was extremely strong evidence to suggest that myofibrils had larger areas and diameters than mitochondria. The myofibrils were wider on average by 0.35-0.42 µm and larger by 0.16-0.20 µm^{2} (p ≈ 0 in both cases). There was also extremely strong evidence that the median diameter varies (p ≈ 0 for both myofibrils and mitochondria).

3.2 Centroid patterns

3.2.1 Organelle pattern characteristics

The fitted intensity models were able to detect the mitochondrion cluster without manual intervention. The presence of a cluster means that a higher order polynomial (cubic or quartic) has a lower AIC value than the homogeneous intensity model. For all of the myofibril and mitochondrion patterns except one, the homogeneous model had a lower AIC value than the cubic and quartic models (for example Figure 7), indicating that the homogeneous model may be sufficient to describe the distribution of centroids within the cell.

(a)

(b)

Figure 7: Intensity models for a homogeneous myofibril pattern, with the corresponding AIC values;

(a) homogeneous intensity (lowest AIC value) and (b) 3rd order polynomial.

However one cell had a large clump of mitochondria - similar to the clusters shown to form at the nuclear poles (??). For this cell’s mitochondrion centroid pattern the cubic model had a lower AIC value (AIC = 60.9) than the homogeneous model (AIC = 68.2), indicating that a higher-order polynomial model is more appropriate than a homogeneous model.

Figure 8: Mitochondrion pattern for which the cubic model (shown) had the lowest AIC value.

The G-function and PCF showed similar interaction between myofibrils and mitochondria in all the cells (Figure 9). All of the patterns display evidence of a hardcore process (no points within a certain r_{min} of each other), and in general there are more myofibril-mitochondrion pairs (green) close together than there are myofibril-myofibril pairs (red) or mitochondrion-mitochondrion pairs (blue). This is evident in Figure 9, where the green lines lie above the red and blue lines at around r ≤ 0.5. G-function and PCF.

Monte Carlo tests for each pattern indicated that the Poisson model should be rejected, as the observed curves lay below the lower envelope.

(a)(b)

Figure 9: Interaction statistics for all the centroid patterns; (a) the nearest neighbour distance cumulative

distribution function (G-function) and (b) the pair correlation function (PCF). Red = distances between myofibrils. Blue = distances between mitochondria. Red = distances from myofibrils to mitochondria (cross function). Black = theoretical curves for a Poisson process (no interaction): G(r) = 1 − e^{−λπr^{2}}

and g(r) = 1

3.2.2 Model results

When the fitted thinned Poisson models were simulated and compared to the orginal data, it became evident that there was a region of repusion beyond the hardcore radius.

Figure 10b shows one simulation of the thinned Poisson model fitted to Figure 10a. The Monte Carlo test (Figure 10c) shows that the observed curve lies below the lower simulation envelope around r = 0.5, so the thinned Poisson model was rejected and a region of interaction beyond the hardcore radius was introduced with the Strausshardcore model. Similar results were obtained for each of they myofibril and mitochondrion patterns. This further repulsion could be due to the fact that the organelles are different sizes. The observed hardcore characteristics are due the fact that the organelles have finite radii, so the minimum distance would be connected to the size of the smallest organelles. However there were also many larger organelles, whose centroids would be unable lie within r_{min} of each other, contributing to an observed region of ’repulsion’ outside the hardcore radius.

The Strauss-hardcore model was able to reproduce the characteristics seen in the centroid pattern (Figure 10e). A visual comparison of the simulated patterns (Figures 10b and 10d) illustrates how the Strauss-hardcore model produces more regular patterns than the thinned Poisson model – similar to the original pattern, which is quite regular (Figure 10a).

(a)

(b) (c)

(d) (e)

Figure 10: Representative results for interaction models fitted to the myofibrils of one cell: (a) original centroid pattern, (b) a simulation of a thinned Poisson model, (c) Monte Carlo test for that model, (d) a simulation of the Strauss-hardcore Gibbs model and (e) Monte Carlo test for that model. In (c) and (e), ‘obs’ = PCF for the observed centroid pattern (Figure 10a), ‘mmean’ = mean curve from the 99 model simulations. ‘hi’ and ‘lo’ =upper and lower boundaries of the simulation envelopes.

Strauss radii were fitted for both myofibrils and mitochondria as single-type patterns, as outlined in Section 2.2.3. The mean Strauss radius for myofibril centroids was 0.55 µm (SD = 0.10 µm) and for mitochondrion centroids it was 0.48 µm (SD = 0.13 µm). There was weak evidence to suggest that the fitted Strauss radius for myofibrils is greater than the mitochondrion Strauss radius for the same cell (p = 0.0715), with an average difference of up to 0.15 µm. It was thought that this could be due to the fact that myofibrils are larger on average than mitochondria, so their centroids cannot come so close to each other. There did not seem to be any clear link between the optimal Strauss radius and the organelle size or minimum distance between centroids.

For the multi-Strauss-hardcore model, three Strauss radii were needed: S_{myo-myo}, S_{mit-mit} and S_{myo-mit}. Initially S_{myo-myo} and S_{mit-mit} were taken to be the same as the optimal values for the singletype patterns, however when these values were used

Figure 11: Myofibril and mitochondrion PCFs with Monte Carlo envelopes calculated from 99 simulations of the multi-Strauss-hardcore model fitted to Cell 1. From left to right, top to bottom: g_{myo}(r), g_{mit}(r), g_{myo,mit}(r),

and g_{myo,•}(r) − g(r). The curve ‘obs’ is the observed PCF, calculated from the Cell 1 organelle data, and ‘mmean’ is the mean curve from the 99 simulations of the multi-Strauss-hardcore model. ‘hi’ and ‘lo’ refer to

the upper and lower boundaries of the model simulation envelopes.

to optimise for S_{myo-mit}, no value could be found that resulted in an acceptable model – that is, one for which the myofibril, mitochondrion and cross Gfunction and PCF simulation envelopes each fully enclosed the corresponding observed function. This meant that the three Strauss radii had to be optimised simultaneously – a computationally expensive process.

As a preliminary exercise, the simultaneous optimisation was carried out for only one cell, to determine whether it could yield the desired results. Acceptable (although not yet optimal) values were found for this cell, and Monte Carlo tests from the resulting model are shown in Figure 11. The Strauss radii for this model were: between myofibrils, S_{myo-myo} = 0.61 µm, between mitochondria, S_{mit-mit} = 0.65 µm and between the two types, S_{myo-mit} = 0.48 µm.

The myofibril Strauss interaction radius was the same as the optimal Strauss radius found for the myofibril-only Strauss-Hardcore model (0.61 µm), but the mitochondrion Strauss interaction radius was greater than the mitochondrion-only value (0.55 µm).

It is also interesting to note that the radius of repulsion is smaller between organelles of different types (0.48 µm) than it is between organelles of the same type (0.61 µm and 0.65 µm), meaning myofibrils are more likely to have mitochondria close to them than myofibrils (and vice versa). This indicates that the Strauss interaction radius may not be simply due to organelle size, but rather due to the fact that myofibrils and mitochondria tend to alternate within the cell. This regularity can be confirmed by visually examining the cell images.

4 Conclusions and Future Work

Electron micrographs were obtained of rat ventricular myocytes in cross section, and used to obtain the positions of each myofibril and mitochondrion centroid. Spatial point pattern statistics were used to analyse the arrangement of the organelles and several spatial point pattern models were investigated.

It was determined that orgnelle centroid patterns can be was approximately homogeneous, except for one cell which had one large clump of mitochondria. Single and cross G-function and PCF are were used to characterise the interaction between organelle centroids – both interactions between organelles of the same type, and organelles of different types. These statistics showed that organelle patterns can be modelled as hardcore process. Simulations of simple hardcore models also indicated a region further repulsion beyond the hardcore radius. Multi-Strauss-hardcore Gibbs models were able to capture the characteristics of the arrangements observed in the cell images.

Organelles of different types tended to be found closer together than organelles of the same type – perhaps due to the alternating pattern observable in the images. Mitochondria tend to be closer together than myofibrils.

In order to make statistical inferences about the cardiac cell structure, the analysis outlined here must performed on a larger sample of cells. Similar studies could be carried out on the arrangement of myofibrils and mitochondria in diseased cells, providing the basis for a quantitative analysis of pathological cell structure.

A complete, three-dimensional cell model must include myofibril and mitochondrion shape information, and also structural information in the longitudinal direction. Once the three-dimensional arrangement and structure of myofibrils and mitochondria has been established, the arrangement of other organelles involved in excitation-contraction coupling (such as the TATS and SR) can be modelled around them to create a complete structural model of the cardiac myocyte. This can then be used as a realistic spatial domain, over which mathematical models of cardiac cell function can be solved. Such computational models could be used to quantify the effect that pathological cell structure has on cell function, offering fresh insight into the mechanisms of heart disease, potentially providing targets for treatment.