\documentclass[preprint,authoryear,12pt]{elsarticle}
\usepackage[latin1]{inputenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}

\usepackage{shortcuts}

\begin{document}
\begin{frontmatter}
\title{Quantifying and Modelling the Spatial Organisation of Myofibrils and Mitochondria in Cardiac Myocytes}
\maketitle
\author{Elizabeth Theakston, Cameron G. Walker, Michael O'Sullivan, Vijay Rajagopal\corref{*}}
%\ead{email address}
%% \ead[url]{home page}
%% \fntext[label2]{}
\cortext[*]{v.rajagopal@auckland.ac.nz}
%% \address{Address\fnref{label3}}
%% \fntext[label3]{}
\begin{abstract}

\end{abstract}

\begin{keyword}
electron microscopy, cardiac cell ultrastructure, spatial statistics, spatial organisation
\end{keyword}
\end{frontmatter}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{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 \citep{Schaper91,medeiros91,chen09}, although giant mitochondria can also form \citep{Kanzaki10}. Myofibrils can become disarranged and degenerated in diabetes \citep{Saito84}, and their function is impaired due to the proliferation of mitochondria . The nuclei can also grow and change shape \citep{Schaper91}. Myocardial infarction causes the TATS to become disorganised and disrupted \citep{quinn03,louch06}, while in dilated cardiomyopathy it has been shown to proliferate \citep{Schaper91}. 

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 \citep{modelling03}, 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 \secref{res-sizedata}, but little has yet been done to quantify the \textit{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 cross-section. We used spatial point patterns to analyse and simulate the arrangement of organelle centroids.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Methods}
\subsection{Tissue processing and data acquisition}
\begin{small}

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\subs{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\subs{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 post-stained 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 $\times$ 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 \textbf{spatial point patterns.}

\end{small}


\subsection{Spatial point pattern statistics and models}
A spatial point process is a random pattern of points in \textit{d}-dimensional space [ref?]. These points can have information associated with them -- called a \textbf{mark}. Marks can either be \textbf{numeric} (for example the radius or area of object the point represents) or \textbf{categorical} (to differentiate between objects of different types). Patterns with categorical marks are called \textbf{multi-type point patterns}  \figrefp{marked-example}. 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 \textbf{intensity} of points (how many points there are within a given area) and by the \textbf{interaction} between points (how the positioning of one point affects the positioning of other points around it) \citep{illian08, Baddeley2008}. 

This section gives a brief  description of the functional aspects of spatial point pattern statistics and modelling techniques. \cite{illian08}, \cite{Eckel08} and \cite{Baddeley2008} 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.  

\begin{figure}[htb]
 \begin{center}
 \includegraphics[width=0.25\textwidth]{./img/marked.pdf}
 \end{center}
 \caption[Multi-Type Point Pattern]{An example of a multi-type point pattern, with categorical marks of two types.}
 \label{fig:marked-example}
\end{figure}




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 \textbf{spatstat} library \citep{baddeley05}.

\subsubsection{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 \figref{dens-hom} 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 \textbf{intensity function}) -- for example the pattern in \figref{dens-quad}, in which the intensity function is a quadratic function of the \textit{x} coordinate.

\begin{figure}[htb]
  \centering
   \subfloat[]{\label{fig:dens-hom}\includegraphics[width=0.22\textwidth]{./img/dens-hom}}\hspace{2mm}
  \subfloat[]{\label{fig:dens-quad}\includegraphics[width=0.22\textwidth]{./img/dens-quad}}\\
  \subfloat[]{\label{fig:constant}\includegraphics[width=0.22\textwidth]{./img/constant}}\hspace{2mm}
  \subfloat[]{\label{fig:quadratic}\includegraphics[width=0.22\textwidth]{./img/quadratic}}
 \caption[Intensity Plots]{Example patterns:
  (a) a pattern generated according to a homogeneous intensity distribution and
  (b) a pattern with an intensity function that is quadratic in \var{x}). Scale bars indicate the intensity in points per unit\sups{2}; windows are both size 1~unit~x~1~unit. (c) and (d) represent the intensity functions in \textit{x} for patterns (a) and (b) respectively.}
  \label{fig:intensity-plots}
\end{figure}

\begin{figure*}[htb]
  \centering
  \subfloat[]{\label{fig:clusterh}\includegraphics[width=0.23\textwidth]{./img/clusterh}}\hspace{5mm}
  \subfloat[]{\label{fig:regularh}\includegraphics[width=0.23\textwidth]{./img/regularh}}\hspace{5mm}
  \subfloat[]{\label{fig:poissonh}\includegraphics[width=0.23\textwidth]{./img/poissonh}}
  \caption{Point patterns with an overall \textit{homogeneous intensity}, illustrating the different types of interaction between points: 
  (a) attraction (clustering),
  (b) repulsion (regularity) and
  (c) no interaction (`Poisson').}
  \label{fig:examplepps}
\end{figure*}

The \textbf{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 \textbf{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 \figref{dens-quad} is shown in \figref{modelling-intensity}. The results for models fitted to the cell data are given in \secref{res-spatstat}.

\begin{figure}[htb]
 \begin{center}
 \includegraphics[width=0.3\textwidth]{./img/intensity-modelling}
 \end{center}
 \caption[Fitting Intensity Models Example]{A polynomial (quadratic) intensity models fitted to the pattern in \figref{dens-quad}.   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.}
 \label{fig:modelling-intensity}
\end{figure}


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



The interaction can be investigated using the  cumulative distribution function for the \textbf{nearest neighbour distances} ( \ref{eq:G}) . This is  called the \textbf{\funcs{G}} and denoted $G(r)$. The nearest neighbour distances are the Euclidean distances between each point and its nearest neighbour. 

\begin{equation} \label{eq:G}
	G(r) = \dfrac{1}{N} \sum_{i=1}^{N} \textbf{1} \left( d(x_{i})\leq r \right) \quad \text{for } r \geq 0
\end{equation}

\noindent where $ d(x_{i}) $ is the distance of point $ x_{i} $ to its nearest neighbour, and  \textbf{1}(\textit{expression}) is the \textbf{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 \textit{r} values, it indicates clustering (more points closer together) and  when it lies below the expected curve it indicates regularity (fewer points close together). \figref{Gfunctions} shows the \var{G}-functions for the example patterns shown in Figures \ref{fig:clusterh}, \ref{fig:poissonh} and \ref{fig:regularh}. \figref{regularG} shows that the pattern in \figref{regularh} is likely to come from a \textbf{hardcore process}, as there are no points within $r = 0.18$ of each other. 

\begin{figure*}[ht]
  \centering
  \subfloat[]{\label{fig:clusterG}\includegraphics[width=0.3\textwidth]{./img/clusterG}}\hspace{2mm}
  \subfloat[]{\label{fig:poissonG}\includegraphics[width=0.3\textwidth]{./img/poissonG}}\hspace{2mm}
  \subfloat[]{\label{fig:regularG}\includegraphics[width=0.3\textwidth]{./img/regularG}}\\
    \subfloat[]{\label{fig:clusterP}\includegraphics[width=0.3\textwidth]{./img/clusterP}}\hspace{2mm}
  \subfloat[]{\label{fig:poissonP}\includegraphics[width=0.3\textwidth]{./img/poissonP}}\hspace{2mm}
  \subfloat[]{\label{fig:regularP}\includegraphics[width=0.3\textwidth]{./img/regularP}}\\

  \caption{\funcs{G} and PCF for the example patterns previously shown in Figures \ref{fig:clusterh}, \ref{fig:poissonh} and \ref{fig:regularh}, showing respectively
  (a) \& (d) clustering (obs above theo),
  (b) \& (e) no interaction (obs close to theo) and
  (c) \& (f) regularity (obs = 0 until $r\approx0.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 \func{G}, $G_{theo}(r) = 1 - e^{-\lambda \pi r^2}$ and for the PCF, $g_{theo}(r) = 1$).}
  \label{fig:Gfunctions}
\end{figure*}

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

\begin{equation} \label{eq:Kest}
	K(r) = \frac{1}{\lambda N} \sum^{N}_{i} \sum^{N}_{j, j\neq i} \mathbf{1}(d_{ij} \leq r)
\end{equation}

\noindent where $ \hat{ \lambda} $ is the estimated intensity, $d_{ij}$ is the distance between points \textit{i} and \textit{j}. For a homogeneous Poisson distribution, the theoretical value for the \textit{K}-function is  $\pi r^2$.

The \textbf{pair correlation function} (PCF), denoted $g(r)$, is the derivative of the \textit{K}-function normalised by $2 \pi r$ (expected curve from a Poisson process).

\begin{equation}\label{eq:pcf}
g(r) = \dfrac{1}{2\pi r} \dif{K}{r}
\end{equation}

The PCF indicates how many pairs of points lie at a distance \var{r} apart, proportional to the expected number for a Poisson process. Thus when \var{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). \figref{Gfunctions} shows the PCFs for the patterns in Figures \ref{fig:clusterh}, \ref{fig:poissonh} and \ref{fig:regularh}. 

%\begin{figure}[ht]
%  \centering
%  \subfloat[]{\label{fig:clusterP}\includegraphics[width=0.2\textwidth]{./img/clusterP}}\hspace{2mm}
%  \subfloat[]{\label{fig:poissonP}\includegraphics[width=0.2\textwidth]{./img/poissonP}}\hspace{2mm}
%  \subfloat[]{\label{fig:regularP}\includegraphics[width=0.2\textwidth]{./img/regularP}}\\
%  \caption[Example Pair Correlation Functions]{Pair correlation functions (PCFs) for the example patterns shown previously in Figures \ref{fig:clusterh}, \ref{fig:poissonh} and \ref{fig:regularh}, showing
%  (a) clustering,
%  (b) no interaction and
%  (c) regularity (graph is zero until $r\approx0.18$, so no pairs of points within that distance of each other) respectively. The curve `obs' is the observed function for the pattern, and the dotted line `theo' is the theoretical Poisson curve ($g(r) = 1$).}
%  \label{fig:PCFs}
%\end{figure}


For \textbf{multi-type marked point patterns}, the \funcs{G} and PCF can be calculated seperately for distances between points of the same type, and also as \textbf{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 (\var{i}) and the point of the other type (\var{j}) closest to it (or vice versa). $K_{ij}(r)$ is the average number of type \var{j} points within \var{r} of a  point of type \var{i}, normalised by the local intensity of type \var{j} points($\lambda_{j}$) -- or vice versa. The cross-PCF, $g_{ij}(r)$ is calculated from $K_{ij}(r)$, as described in (\ref{eq:pcf}).

In this study, the \funcs{G} 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 \figref{interaction} in \secref{res-spatstat}.

Due to the stochastic nature of point processes, it is helpful to compare observed functions to a \textit{range} of acceptable values rather than a single theoretical curve. This can be done using \textbf{Monte Carlo envelopes} -- found by generating 99 simulations of a homogeneous Poisson process, calculating the corresponding \funcs{G} or PCF curves, and taking the highest and lowest  function values for each value of \textit{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 an observed pattern is consistent with a particular model, by calculating Monte Carlo envelopes from 99 simulations of that model.

\subsubsection{Modelling interaction}\label{sec:modelling}
Organelles have a finite radius, so the centroids cannot come within a certain distance \rmins of each other. This means they must be modelled with  a \textbf{hardcore process}. A simple hard-core model is the \textbf{thinned Poisson process}. Points are \textit{generated} according to the Poisson process (no interaction), but if they fall within  the hardcore radius (\rmin) of another point, they are \textit{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 \funcs{G} and PCF. 

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


\[
\phi(r) = \left\{ 
\begin{array}{l l}
  \infty & \quad r \leq r_{hc}\\
  \beta & \quad r_{hc} < r \leq r_{S}\\
  0 & \quad r > r_{S}
\end{array} \right.
\]

\begin{figure}[htb]
 \begin{center}
 \includegraphics[width=0.23\textwidth]{./img/PPF}
 \end{center}
 \caption{Pair potential function for a Strauss-hardcore model.}
 \label{fig:ppf}
\end{figure}

A plot of $ \phi (r) $ for the Strauss-hardcore process is shown in \figref{ppf}. This pair potential function means that below a certain radius there is \textit{infinite repulsion}, so a point cannot fall within a hardcore radius $ r_{hc} $ of another.  The range  $r_{hc} < r \leq r_{S} $ defines the region of \textbf{Strauss interaction}. Points within this distance of each other will experience repulsion, the strength of which is defined by $\beta$.  In theory,  $\beta$ 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 \cite{Baddeley2008} and Section 6.6 in \cite{illian08}. 

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 ($\beta$) 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 \times$ \rmins (\rmins =  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 \textbf{multi-Strauss-hardcore} process. The latter required finding separate Strauss radius values for myofibril-myofibril, mitochondrion-mitochondrion and myofibril-mitochondrion interactions. Hardcore radii were set to 0.98 $\times$ \rmin.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{table*}[htb] \caption[Cell and Organelle Size Statistics]{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.} \label{tab:cellstats}
	\center\small
	\begin{minipage}{1\textwidth}
	\center\small
	\renewcommand{\thefootnote}{\arabic{footnote}}
	\begin{tabular}{l p{1.5cm} l p{1.6cm}p{1.5cm} cr@{$\pm$}l cr@{$\pm$}l}
	\toprule
	\multirow{2}{1.5cm}{Cell area (\textmu m)}   
	 	& \multirow{2}{1.5cm}{Organelle} 
		 	&  \multirow{2}{1.5cm}{Density (\pmicronsq)}
		 		& \multirow{2}{1.7cm}{\% of Cell area}
	 				& \multicolumn{3}{c}{Diameter (\micron)}
	 					& \multicolumn{3}{c}{Area (\micronsq)} \\
	%	\cmidrule(lr){7-9}\cmidrule(lr){10-12} 
		& & & & Med & Mean & \%SD & Med & Mean & \%SD \\
	\midrule
 	 \multirow{4}{*}{63.5$\pm$45\%} 
				& \multirow{2}{*}{myofibrils} & 1.42$\pm$19\% & 50.4 $\pm$7\%& 1.02 & 1.11 & 43\% & 0.30 & 0.37 & 69\%\\
				& & \grey{-}
						& \grey{54-61}\footnote{\cite{anversa85,barth92,crisman85,frenzel88,medeiros91,schaper85,skepper95}}
								& \multicolumn{2}{c}{\grey{0.7-1.1}\footnote{\cite{hansen78,linke94,stehle02}}} 
										& & \multicolumn{2}{c}{\grey{-}}\\
 	   	   	   	& \multirow{2}{*}{mitochondria} &  1.92 $\pm$20\%& 33.9$\pm$10\% & 0.68 & 0.73 & 48\% & 0.15 & 0.19 & 81\%\\
				& & \grey{approx. 1}\footnote{\cite{chen09}}
						& \grey{26-36}\footnote{\cite{anversa85,barth92,crisman85,frenzel88,medeiros91,oron85,schaper85,skepper95}}
								& \multicolumn{2}{c}{\grey{0.5-1}\footnote{\cite{scheffler99}}}
										& & \multicolumn{2}{c}{\grey{0.37,0.57}\footnote{\cite{hom10} and \cite{chen09} respectively}}\\
	\bottomrule
	\end{tabular} \par
   \vspace{-0.75\skip\footins}
   \renewcommand{\footnoterule}{}
  \end{minipage}
\end{table*}



\section{Results and Discussion}


\subsection{Size data}\label{sec:res-sizedata}
The average statistics for cell and organelle sizes are shown in \tabref{cellstats}, 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\microns and larger by 0.16-0.20\micronsqs ($p\approx 0$ in both cases). There was also extremely strong evidence that the median diameter varies  ($p\approx 0$ for both myofibrils and mitochondria).  

\subsection{Centroid patterns}

\subsubsection{Organelle pattern characteristics}\label{sec:res-spatstat}

The fitted \textbf{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 \figref{intensitymodels}), indicating that the homogeneous model may be sufficient to describe the distribution of centroids within the cell.

\begin{figure}[htb]
  \centering
  \subfloat[]{\label{fig:cell1mfs0}\includegraphics[width=0.30\textwidth]{./img/cell1mfs0}}\\
  \subfloat[]{\label{fig:cell1mfs3}\includegraphics[width=0.32\textwidth]{./img/cell1mfs3}}
  \caption[Intensity models]{Intensity models for a homogeneous myofibril pattern, with the corresponding AIC values; 
  (a) homogeneous intensity (lowest AIC value) and
  (b) 3rd order polynomial.}
  \label{fig:intensitymodels}
\end{figure}

However one cell had a large clump of mitochondria - similar to the clusters shown to form at the nuclear poles \citep{scheuermann93,hom09}. 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. 

\begin{figure}[htb]
 \begin{center}
 \includegraphics[width=0.45\textwidth]{./img/cell3mit3}
 \end{center}
 \caption{Mitochondrion pattern for which the cubic model (shown) had the lowest AIC value.}
 \label{fig:cell3mit3}
\end{figure}


The \textbf{\funcs{G}} and \textbf{PCF} showed similar interaction between myofibrils and mitochondria in all the cells \figrefp{interaction}. All of the patterns display evidence of a hardcore process (no points within a certain \rmins 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 \figref{interaction}, where the green lines lie above the red and blue lines at around $r\leq 0.5$. \funcs{G} and PCF.


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

\begin{figure}[H]
  \centering
  \subfloat[]{\label{fig:G}\includegraphics[width=0.295\textwidth]{./img/Gest}}\hspace{2mm}
  \subfloat[]{\label{fig:P}\includegraphics[width=0.295\textwidth]{./img/PCF}}
  \caption[Interaction statistics]{Interaction statistics for all the centroid patterns; 
  (a) the nearest neighbour distance cumulative distribution function (\func{G}) 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^{-\lambda \pi r^2}$ and $g(r) = 1$}
  \label{fig:interaction}
\end{figure}

\subsubsection{Model results}

When the fitted \textbf{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. 


\figref{SSI-sim} shows one simulation of the thinned Poisson model fitted to \figref{Cell1mfs}. The Monte Carlo test \figrefp{SSI-pcf} shows that the observed curve lies below the lower simulation envelope around \var{r} = 0.5, so the thinned Poisson model was rejected and a region of interaction beyond the hardcore radius was introduced with the Strauss-hardcore 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 \rmins of each other, contributing to an observed region of 'repulsion' outside the hardcore radius.

The \textbf{Strauss-hardcore} model was able to reproduce the characteristics seen in the centroid pattern \figrefp{SH-pcf}. A visual comparison of the simulated patterns (Figures~\ref{fig:SSI-sim} and \ref{fig:SH-sim}) illustrates how the Strauss-hardcore model produces more regular patterns than the thinned Poisson model -- similar to the original pattern, which is quite regular \figrefp{Cell1mfs}.

\begin{figure}[htb]
  \centering
  \subfloat[]{\label{fig:Cell1mfs}\includegraphics[width=0.21\textwidth]{./img/Cell1mfs}}\\
  \subfloat[]{\label{fig:SSI-sim}\includegraphics[width=0.225\textwidth]{./img/SSI-sim}}\hspace{1mm}
  \subfloat[]{\label{fig:SSI-pcf}\includegraphics[width=0.225\textwidth]{./img/SSI-pcf}}\\
  \subfloat[]{\label{fig:SH-sim}\includegraphics[width=0.225\textwidth]{./img/SH-sim}}\hspace{1mm}
  \subfloat[]{\label{fig:SH-pcf}\includegraphics[width=0.225\textwidth]{./img/SH-pcf}}
  \caption[Thinned Poisson and Strauss-Hardcore Gibbs Simulation Results]{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 \figrefp{Cell1mfs},  `mmean' = mean curve from the 99 model simulations. `hi' and `lo' =upper and lower boundaries of the simulation envelopes.}
  \label{fig:modelling-single}
\end{figure}

Strauss radii were fitted for both myofibrils and mitochondria as \textbf{single-type} patterns, as outlined in \secref{modelling}. The mean Strauss radius for myofibril centroids was  0.55\microns (SD = 0.10\micron) and for mitochondrion centroids it was 0.48\microns (SD = 0.13\micron). There was weak evidence to suggest that the fitted Strauss radius for myofibrils is greater than  the mitochondrion Strauss radius for the same cell (\textit{p} = 0.0715), with an average difference of up to 0.15\micron. 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.

\begin{figure*}[ht]
 \begin{center}
 \includegraphics[width=0.85\textwidth]{./img/cell1_MSH}
 \end{center}
 \caption[Multi-Strauss-Hardcore Simulation Results]{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,\bullet}(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.}
 \label{fig:res-MSH}
\end{figure*}

For the \textbf{multi-Strauss-hardcore} model, three Strauss radii were needed: $ S_{\text{myo-myo}}$, $ S_{\text{mit-mit}} $ and $ S_{\text{myo-mit}} $. Initially $ S_{\text{myo-myo}}$ and $ S_{\text{mit-mit}} $ were taken to be the same as the optimal values for the single-type patterns, however when these values were used to optimise for $ S_{\text{myo-mit}} $, no value could be found that resulted in an acceptable model -- that is, one for which the myofibril, mitochondrion and cross \funcs{G} 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 \figref{res-MSH}. The Strauss radii for this model were: between myofibrils, $S_{\text{myo-myo}} = 0.61$\micron,  between mitochondria, $S_{\text{mit-mit}}=0.65$\microns and between the two types, $S_{\text{myo-mit}} =0.48$\micron. 

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

It is also interesting to note that the radius of repulsion is smaller between organelles of different types (0.48\micron) than it is between organelles of the same type (0.61\microns and 0.65\micron), 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.
 



\section{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.



\renewcommand{\bibname}{References}
\footnotesize
\bibliographystyle{mybibstyle}
\bibpunct{(}{)}{,}{a}{,}{,}
\bibliography{bibliography.bib}


\end{document}
