ImageAnalStereol2009;28:11-22 OriginalResearchPaper A STEREOLOGICALUNFOLDINGMETHODFOR THESTUDY OF THE MITOCHONDRIALNETWORK DANIELCHRISTOPHER FLYNN1,JAMESDARWINNULTON2,PETERSALAMON2, TERRENCE GENEFREY3,AVINOAM RABINOVITCH4,MEIGUOSUN3ANDARLETTEROSE CORNELIA BALJON1 1DepartmentofPhysics,SanDiegoState University,SanDieg o,California,USA;2DepartmentofMathematical Sciences,SanDiegoState University;3DepartmentofBiology,SanDiegoState University;4Departmentof Physics,Ben-GurionUniversityofthe Negev,Beer-Sheva,I srael e-mail: mdacflynn@umich.edu,jnulton@mail.sdsu.edu,sal amon@math.sdsu.edu,tfrey@sciences.sdsu.edu, sun.mei@gene.com,abaljon@mail.sdsu.edu (AcceptedMarch6,2009) ABSTRACT When the mitochondrialnetworkwithin the cell is in a fragme ntedstate, it resemblesa spatial distribution of prolatespheroidsofvariousshapes,sizesandorientation s. Thispaperpresentsamaximum-likelihoodscheme for inferring the distribution of spheroidal shapes, sizes and orientations from the observed distribution of elliptical sections. We also present a parametric bootstra p suitable for inferring confidence intervals for the parametersdescribingtheshapeandsize distributionofth espheroids. Keywords:bootstrap,maximumlikelihood,mitochondria,s pheroid,statistics, unfoldingproblem. INTRODUCTION Themotivationforthispapercomesfromthestudy of conformationalchangesin mitochondria associated with cell death. During a late stage of apoptosis (programmed cell death) the cell’s mitochondrial network fragments into individual mitochondria. Cross-sections of these particles observed in electron micrographs suggest that they have an oblong shape. The particles originate in a process that is reminiscent of the way sausage is made, in that the fragments apparently result from the pinching off (fissioning) the cylindrical tubes that make up the mitochondrial network.Sincethetubesarecylindrical,weexpectthe fragments, being sacks of fluid, to adjust to a roughly spheroid shape. In our study we approximate these particles by prolate spheroids. Low resolution visual evidencefromthree-dimensionalconfocalmicroscopy (Franket al., 2001) argues against the presence of oblate structures. Moreover, in section “A sample application to mitochondria” we show that the sectioneddataitself arguessimilarly. These particles are believed to undergo changes over time in shape and size due to disruption of inner membrane function and possible resultant swelling. While confocal light microscopy is able to partially document this process, the resolution is limited to values comparable to the particle size. Quantitativeinformationonsizeandshapeisobtained from electron micrographs, which usually show only 2-dimensional cross-sections of the particles.Nevertheless, assuming the particles are prolate spheroids, it is possible to infer 3-dimensional information from the statistics of the 2-dimensional sections. It is this stereological problem that we propose to address here. Wicksell (1925; 1926) dubbed this “the corpuscleproblem.”He was motivatedby the studyof secondary follicles embedded in the lymphatic tissue of the spleen. More recently the problem has become known as the “unfolding problem.” An excellent history and set of references can be found in Stoyan et al.(1987). A complete treatment of stereological methods with an emphasis on practical problems can be found in the two volumes of Weibel (1979; 1980). Mostrecently,thegeneralproblem,involving3Dsize, shapeandorientation,wastreatedbySato etal.(1996). The basic mathematical problem is to relate the distributions of observable quantities in 2D cross- sections to the distributions of desired quantities that are descriptive of the 3D particles. Typically this takes the form of an integral transform expressing, for example, the distribution Ψof cross-sectional diameters in terms of the distribution Φof particle diameters. The practical problem is then to use this machinery to work from data in the form of a sample fromΨto a description of Φ. Two approaches were proposed in Wicksell’s original papers and have been takenupby subsequentauthors(De Hoff, 1962;Cruz- orive, 1976; Franklin, 1977; Sato et al., 1996). First, theintegraltransformcanbeusedtoderivearecursive 11 FLYNNDCET AL:Stereologicalunfoldingmethod scheme that expresses the moments of Φin terms of those of Ψ. Secondly, in some cases, it is possible to invert analytically the transform and to solve for Φin terms of Ψand its derivatives. Both of these approachessufferwhenthedatasetisnotlargeenough to provide a clear picture of Ψ. In the first approach, statisticalmomentsfromthedataareonlyestimatesof Ψ’s true moments, and the associated statistical error is propagated through the transform. Similarly, in the second approach, estimates of the derivatives of Ψby finite difference methods introduce uncertainties that aredifficultto quantify. Four distinct additional approaches have been identified (Bl¨ odner et al., 1984); among these is the parametric method that we propose to apply. The parametric approach has been used to predict a very precise distribution of sphere diameters from circular cross-sections (Keiding and Jensen, 1972). In the present paper we describe a parametric approach, which does not seem to have been used by previous authors in the context of the full unfolding problem described by Sato, involving simultaneous predictions of the distributions of size, shape and orientation. We propose parametric models for desired distributions forΦ, and use a maximum-likelihood approach to optimize the parameters. The method of parametric bootstrap (Efron and Tibshirani, 1993) was used to determineconfidenceregionsforthe parameters. In the next section we set out the geometrical framework, define the quantities studied, and state the relations among them. In the section following that we apply the method of Sato to derive the machinery expressing the distributions of cross- sectional size and shape in terms of the hypothetical or desired distributions of particle size, shape, and orientation. The three sections beginning with the Maximum-likelihood Scheme present a description of our method. We describe the parametric models used for the hypothetical distributions and the maximum- likelihood scheme used to determine values for the parameters. We also describe how to set up a parametric bootstrap for determining confidence regions for the values of the parameters. The last section presents the results of applying the method to a setofdataderivedfrom electronmicrographs. GEOMETRICALPRELIMINARIES In this section we summarize the geometrical relationsamongtheslicingplane,thespheroid,andthe elliptical cross-section.SeeFig.1.crb a sh 1234 5 θ Fig.1.Drawingofslicedspheroid:1)Axisofspheroid. 2)Prolatespheroid.3)Slicingplane.4)Tangentplane. 5)Elliptical crosssection. We take the slicing plane to be horizontal, a distancesfromthecenterofthespheroid.Wetake θto be the angle between the axis of the spheroid and the vertical. The spheroid is described by its semi-minor axisr, its semi-major axis c, and its eccentricity e0. Similarly, the elliptical section, is described by axes b anda, and eccentricity e. The axes and eccentricities are related by e2 0=1−r2/c2ande2=1−b2/a2. The figure shows a cross-section through the spheroid’s axis. The following equations are easily established and summarize the relationships among the variables requiredforthe analysisin thenextsection. e=e0sinθ, (1) b=r/radicalbigg 1−s2 h2, (2) h=r/radicalBig 1−e2 0sin2θ /radicalBig 1−e2 0. (3) Geometrically θcanbeconsideredtovaryfrom0to π, ands, from −∞to+∞1. Nevertheless,the distribution ofθis understood to be symmetric about the value π/2. Similarly, sis distributed symmetrically about 0. These symmetry considerations allow us, in the following analysis, to restrict θto[0,π/2]andsto [0,∞). RELATIONSAMONGTHE DISTRIBUTIONS Our goal is to extract information on the distributionsof e0,r,andθfromobserveddistributions 1Ofcourse anellipticalslice isobtained only when |s|bande0>e, I(b,e,r,e0) =b√ 1−e2Φ0(sin−1(e/e0)) √ r2−b2/radicalBig 1−e2 0/radicalBig e2 0−e2,(13) andthevaluezeroforallothervaluesof rande0. Thisresult,togetherwiththelimitsonthevariables rande0,canbesubstitutedintoEq.6andtheresulting doubleintegralfactoredto give Ψ(b,e)∝ Ψ1(b)·Ψ2(e), (14) where Ψ1(b) =/integraldisplay∞ bbΦ1(r)√ r2−b2dr, (15) and Ψ2(e) =/integraldisplay1 e√ 1−e2Φ0(sin−1(e/e0))Φ2(e0)/radicalBig 1−e2 0/radicalBig e2 0−e2de0. (16) We caution the reader that Ψ1andΨ2are not probability densities, because they have not been normalized. Normalization is not required to set up the likelihood functions described in the next section. However, as we will see later, these functions are the basic tools for simulating a slicing experiment. In that casetheymustbenormalizedbeforetheycanplaythat role. THE MAXIMUM-LIKELIHOOD SCHEMEANDPARAMETRIC MODELS We have yet to assign statistical models for the densities Φ0(θ),Φ1(r), andΦ2(e0). Once that is done the parameters in those models become the basis for varying the hypothetical distribution of spheroids so as to optimally fit the 2D data. In this section we first formulate the likelihood of a dataset given the hypothetical set of parameter- dependent Φ’s. Maximizing the likelihood optimizes the fit and completes the prediction. We then turn to a description of and rationale for the statistical models recommended for the Φ’s. These were chosen for the applicationdescribedinthe lastsection. Suppose we have collected data in the form of melliptical cross-sections. For each section we haverecorded the semi-minor axis band the eccentricity e:(b1,e1),(b2,e2),...,(bm,em). The data is separated intoD1={b1,b2,...,bm}andD2={e1,e2,...,em}2. D1will be used to infer Φ1(r)through Eq. 15, and D2will be used separately to infer Φ0(θ)andΦ2(e0) throughEq.16. If the distributions Φ0,Φ1, andΦ2in Eqs. 15 and 16 were a true description of the population of spheroids, then D1would be a collection of m independentchoices from the distribution Ψ1(b), and, similarly, D2, a collection from Ψ2(e). Therefore, we can formulate the likelihood of each data set as a product: L1(D1) =m ∏ j=1Ψ1(bj),L2(D2) =m ∏ j=1Ψ2(ej). (17) Maximizing these likelihoods over the space of parameters on which the densities Φ0(θ),Φ1(r), and Φ2(e0)depend determines values of those parameters thatoptimallyfittheobserveddata.Thiscompletesthe algorithm by which predictions for the distribution of spheroids is made. For example, Φ1will be modelled by a normal distribution with the usual parameters µandσ. Then, for the given data D1,L1becomes a function of µandσ. The values µ∗andσ∗that maximize L1then describe the most likely normal model for Φ1that is consistent with the data D1. A similarprocesswith L2isusedtoobtainpredictionsof Φ0(θ)andΦ2(e0)from the data D2. We now describe and discuss the parameter- dependent statistical models recommended for the Φ’s. In each case the simplest appropriate model was chosen. We do not regard these models as strictly valid,orashavingbeenjustifiedbysomeindependent experiment.Rather, they shouldbe regardedas probes of the first- and second-order features of the statistics wearetrying to study. The statistical model suitable for probing the first two moments of an unknown statistics is the normal or Gaussian density; we find that appropriate for Φ1(r). The usual 2-parameter choice for modelling the statistics of a naturally bounded variable is the beta distribution or density; we find that appropriate for bothΦ0(θ)andΦ2(e0), since theta is confined to the interval [0,π]ande0is confined to the interval [0, 1].Φ0requires further comments, given below. There is a symmetry that sets the mean of θatπ/2. As a result Φ0becomes a 1-parameter density. The remaining parameter characterizes the spread, which, asexplainedbelow, canbeinterpreted as a measureof thedeviationfrom isotropic orientation. 2As we have remarked, for the data we collected, there was no di scernable correlation between bande. That result is consistent with equation (14), and, thus, withour assumption that rande0for the spheroids are independent. Hence, the separation of the data. 14 ImageAnalStereol2009;28:11-22 Concretely, the models that we propose are as follows.Φ1is modeledby anormaldensity. Φ1(r|µ,σ) =1√ 2πσexp/parenleftbigg −1 2(r−µ)2 σ2/parenrightbigg .(18) Φ2is modelledby abetadensity, Φ2(e0|p,q) =Γ(p+q) Γ(p)Γ(q)ep−1 0(1−e0)q−1,(19) for whichthemeanandvarianceare givenby µ=p p+q,σ2=pq (p+q)2(p+q+1).(20) For the density Φ0, describing the orientations of the spheroids, we use a simple 1-parameter beta density scaledtothe interval [0,π]: Φ0(θ|α) =2Γ(2α+2) Γ(α+1)2π2α+1θα(π−θ)α.(21) The case of θ, the angle between the spheroid’s major axis and the vertical, requires more discussion because of the special role it plays in the type of data we are studying. First of all, the mitochondrial spheroids are presumably floating freely in the cytoplasm, so why assume anything other than a purely isotropic orientation for their major axes? The micrographsarehorizontalsectionsthroughacellthat isadsorbedtoahorizontalplanesurface.Thereforewe have no reason to expect anisotropy in the horizontal plane.However,thereisreasontoexpectanisotropyin the vertical orientation. The basis of this expectation is related to the origin of the spheroids. The spheroids result from the fragmentation of an extensive tubular networkasexplainedintheintroduction.Sincethecell is strongly attached to a horizontal surface, there is a flattening of the cell from its normal state. Confocal microscopy, our own as well as that of Frank et al. (2001), demonstrate the flattening and indicate that the tubes in the network show some preference for the horizontal. As a result the tubes in the network may show some preference for the horizontal. When the tubes fission, producing the spheroids, a trace of that preference may linger in the orientation of the spheroids. For that reason we structure our model so that it has the potential to detect and quantify any anisotropythatmightexist. The ability to quantify this anisotropy is not the most important reason for including the distribution of spheroid orientation in our scheme. If anisotropy is present and we do not account for it in our likelihood machinery, our results on the distribution of spheroid shape will be distorted. Eq. 1 reveals the intimateconnection between the shape of the elliptical section and the shape and orientation of the spheroid. Eq. 16 shows how their distributions are intertwined in the unfoldingprocessfrom 2D to3D. We now explain our choice of Φ0and how we propose to quantify the anisotropy. For a purely isotropic distribution of points on a sphere, the colatitudes(correspondingtotheangle θ)ofthepoints havethedensityfunction Φ(θ)=1 2sinθon[0,π].This densityhasaneasilycalculatedvarianceof (π2−8)/4. Avariancesmallerthanthisisindicativeofsomething called a “girdle” or “equatorial” distribution. As presentedintheliterature(Selby,1964;Watson,1965), itisadistributionofpointswithahigherconcentration near the equator, and one that is symmetric about the equator. This is the situation we find when considering the orientation of spheroid particles with a slight preference for the horizontal. In our case, the symmetry is a built-in artifice, i.e., a spheroid whose axis is not perfectly horizontal is represented both by the angle θandπ−θ. Weexpectourmodel Φ0toquantifytheanisotropy in the following way. Let σ2be the variance of the particular Φ0determined by our maximum-likelihood algorithm. Wedefineanindexofisotropy: Iiso=σ σiso=2σ√ π2−8. (22) This has a ready interpretation as the ratio of the observed spread in the polar angle to the isotropic spread.Ithasavaluerangingfrom0to1forequatorial distributions. All one needs in order to capture this anisotropy isafairlyall-purpose,1-parameter,symmetricdensity. As reported above our choice was a beta distribution. Other 1-parameter models for equatorial distributions havebeenproposed,andwouldprobablyserveequally well as models for Φ0. For example, one prominent choice among these is the Dimroth-Watson model (Mardia, 1972;Batschelet,1981): Φ(θ) =Ce−κcos2θsinθ, (23) whereCis a suitable normalization constant For κ≥0 it represents equatorial distributions. Here the parameter κcan be regarded as a measure of anisotropy. However, except for κ=0, which corresponds to isotropic, its numerical value has no compelling interpretation. One could claim that this model has the advantage that the isotropic case is included in the family. However, that is not a practical advantageoverEq. 21. Remarkably, for α= 1.16072, Eq. 21 is virtually indistinguishable from sinθ, uniformly over the interval [0,π], and, hence includesthe isotropiccaseaswell. 15 FLYNNDCET AL:Stereologicalunfoldingmethod THE PREDICTORANDTHE SIMULATOR Inthefirstthreeparagraphsoftheprevioussection we have described an algorithm for which the input is data (b1,e1),...,(bm,em)in the form of observed elliptical sizes and shapes, and for which the output is a set of predicted values (µ∗,σ∗,p∗,q∗,α∗)for the parameters that describe the most likely population of spheroids to have produced the input data in a slicing experiment.Letuscallthisalgorithmthe PREDICTOR . As a result of the factorization in Eq. 14, there are actually two algorithms: PREDICTOR 1that predicts the parameters (µ∗,σ∗)for the size distribution of prolate spheroidsfrom the data D1= (b1,...,bm),and PREDICTOR 2that predicts the parameters (p∗,q∗,α∗) for the shape and orientation distributions from the dataD2= (e1,...,em). Notethatthestatisticalroleofthe PREDICTOR isto provide a most likely point estimate for the unknown parameters consistent with a given set of data. A completestatisticalanalysisrequires,inaddition,some measure of the uncertainty of the method in the form of a region, containing the point estimate, within which the true values of the parameters are believed to fall with some specified level of confidence. For that purpose it will be important to define a second algorithm, the SIMULATOR , which uses the same apparatus as the PREDICTOR , but which, in a sense, reversesthe process. The geometry and statistics of the ideal slicing experiment has already been formulated in the mathematical machinery of the two previous sections. Suppose that we know (or hypothesize) that the spheroid population is described by the parameter values (µ0,σ0,p0,q0,α0)3. These values specify the distributions Φ0(θ),Φ1(r),andΦ2(e0)bywayofEqs. 18-21. Eqs. 15 and 16 then uniquely determine the distributions Ψ1(b)andΨ2(e), which describe how we expect any data obtained from the hypothesized population to be distributed. As remarked earlier, Ψ1 andΨ2are presented in non-normalized form, which makes them suitable for a likelihood function, but not suitable as probability densities. Let ˆΨ1(b)and ˆΨ2(e)represent these distributions after they have been normalized so that their integrals are unity. These densities are the tools by which the ideal slicing experiment can be simulated. A simulation producing melliptical sectionsiseffectedbychoosing numbers (b1,...,bm)according to the distribution ˆΨ1(b), and numbers (e1,...,em)from the distributionˆΨ2(e). Thus, the SIMULATOR is the algorithm whose inputis a setof hypotheticalvalues (µ0,σ0,p0,q0,α0) and whose output is a set of data (b1,...,bm) and(e1,...,em)representing the result of an ideal slicing experiment with the population of spheroids whose distribution is described by the parameters (µ0,σ0,p0,q0,α0). Toseehowthe SIMULATOR bearsontheproblems ofreliability andthe accuracyofourmethod,consider the following. Suppose the SIMULATOR is applied to the hypothetical values (µ0,σ0,p0,q0,α0)N times (a large number) in a Monte Carlo scheme. Each of the Nsimulated data sets can be fed to thePREDICTOR to produce Ndifferent predictions (µ1,σ1,p1,q1,α1),...,(µN,σN,pN,qN,αN). If our method is highly accurate we would expect the points (µ1,σ1,p1,q1,α1),...,(µN,σN,pN,qN,αN) to be closely clustered about the point (µ0,σ0,p0, q0,α0)representing the “true” distribution. Indeed, the uncertainty of our method will be measured by quantifyingthis clustering. Wewilladdressthisuncertaintyinthenextsection, but first, let us point out that, just as the PREDICTOR splits into two algorithms, the same is true of the SIMULATOR . That is, SIMULATOR 1takes as input hypothetical values (µ0,σ0)for the size distribution of spheroids, and, in a Monte Carlo scheme, produces a large number Nof new sets of data, which in turn produceNpredictions (µ1,σ1),...,(µN,σN).Thiscan bepicturedasatwo-dimensionalscattergramofpoints whose center of mass should closely approximate the point(µ0,σ0). Fig. 8 shows that scattergram for the examplein the lastsection.Similarly, SIMULATOR 2in a Monte Carlo scheme produces a three-dimensional scatterofpoints (p1,q1,α1),...,(pN,qN,αN)from the hypotheticaldescription (p0,q0,α0)ofthedistribution ofshapesandorientationsofthespheroids. QUANTIFYINGTHE UNCERTAINTY As remarked above, we must establish a region within which we believe the true values of the parameters will fall with some specified level of confidence. The method described here is called a parametric bootstrap (Efron and Tibshirani, 1993); the essential idea is as follows. Numerous samples (of sizem) are simulated from the ideal distribution represented by our point estimates, µ∗,σ∗, etc., of the parameters. Maximum-likelihood estimates for µand 3In the next section we will take these hypothetical values to be the maximum-likelihood values µ∗,σ∗, etc., predicted from the observed data. 16 ImageAnalStereol2009;28:11-22 σ, etc. are then calculated for each of these samples. The statistics of these estimates provide a basis for determining approximate confidence regions for the originalpointestimates. An ordinary variant of bootstrapping may also be used to obtain robust confidence regions. In that case, the numerous samples required are obtained not from the ideal distribution, but by sampling with replacement many times from the original empirical sample. We chose to use the parametric variant outlined aboveanddescribedbelow. Although the method can be applied to any confidence level, for the sake of definiteness, we take thatlevelto be90%confidenceforthis discussionand for theexamplepresentedin thelastsection. Also, for purposes of describing the method we will focus on the size problem, using PREDICTOR 1 andSIMULATOR 1in a Monte Carlo scheme to describe the two-dimensionalconfidenceregion in the space (µ,σ). The case of shape and orientation is conceptuallysimilar,butthediscussionintheexample will reflect the three-dimensional character of the region. Now the only tool we have for assessing a given prediction is the likelihood function L1in Eq. 17 obtained from our real set of data (b1,...,bm).L1is regarded as a function of the parameters (µ,σ)of the prediction. It is convenient to standardize this as a likelihood ratio L1(µ∗,σ∗)/L1(µ,σ),where (µ∗,σ∗) is the maximum-likelihood prediction. As a practical matterintheimplementationofthemethodweusethe logarithm of the likelihood function, so we define the log-likelihood-difference( LLD) by LLD(µ,σ) =ln/parenleftbiggL1(µ∗,σ∗) L1(µ,σ)/parenrightbigg =ln(L1(µ∗,σ∗))−ln(L1(µ,σ)).(24) We can now define the 90% confidence region in the(µ,σ)-plane. Our analysis is of a single data set (b1,...,bm)for which we have a single maximum- likelihood prediction (µ∗,σ∗). We consider that the experiment of collecting this data set is a relatively expensive affair. Nevertheless, suppose, contrary to fact, that we were able to collect many, say N, real data sets, all of the same size m, collectively producingmanypredictions (µ1,σ1),...,(µN,σN).Let xrepresent the 90th percentile of all the numbers LLD(µ1,σ1),...,LLD(µN,σN). Clearly, the inequality LLD(µ,σ)≤xestablishes a region Rof the (µ,σ)- planewiththepropertythat,withprobability0.90,any repetition of our experiment will result in a predictionthatfallswithintheregion R.Thisispreciselywhatwe meanbya90% confidenceregion. Unfortunately, x, as defined above, is not accessiblewithoutrepeatingourexpensiveexperiment Ntimes.Nevertheless,wecansimulatetheexperiment as many times as we like by taking the most likely prediction (µ∗,σ∗)from our single experiment as the inputto SIMULATOR 1.Sinceinthiscasewearetaking onlyourmostlikelyprediction (µ∗,σ∗)asourestimate ofthe“true”values— (µ0,σ0),this will produceonly anapproximationoftheconfidenceregion. To summarize, our approximation to the 90% confidence region is determined as follows. The maximum-likelihoodprediction (µ∗,σ∗)fromthedata is supplied as input to SIMULATOR 1to produce N sets of simulated data. PREDICTOR 1then produces a prediction for each of those sets of data: (µ1,σ1), ...,(µN,σN). The 90% confidence region is then determined by the inequality LLD(µ,σ)≤ x, where xis the 90th percentile of the numbers LLD(µ1,σ1),...,LLD(µN,σN). A SAMPLEAPPLICATIONTO MITOCHONDRIA InthesectionontheMaximum-likelihoodScheme we discussed our reasons for believing that the 3D structures we are studying exhibit anisotropy in vertical orientation. However our evidence for it comes from three-dimensional confocal microscopy and is visual rather that quantitative. It would have been desirable to have sectioning from a number of different angles (besides just horizontal) in order to quantify the anisotropy directly from the sectioned data.However,thecellswhosemitochondriaarebeing sectioned are grown on a flat surface and are very thin. Non-horizontalsectioningin this circumstanceis technically difficult; in fact for this study, prohibited. Nevertheless, the parameter αin our statistical model provides us with an indirect confirmation of the presence of anisotropy as well as a quantitative measure. In the introduction we promised a direct empirical argument in favor of the prolate hypothesis for the shape of the spheroid particles. Weibel (1979) states tworoughempiricalteststhatcanbeperformedonthe sectioned data to decide between prolate and oblate. As a first test, consider the diameter, d, of the profile of largest area among the 5% most circular profiles. Forcomparison,considerthe largestdiameter, 2 a,and the smallest diameter, 2 b, of the profile of largest area among the 5% most elongated profiles. Weibel agrues 17 FLYNNDCET AL:Stereologicalunfoldingmethod thatifd≈2b,the3Dstructuresaremostlikelyprolate; ifd≈2a, they are most likely oblate. In our case, d= 340 nm, 2 a= 890 nm, and 2 b= 270 nm, pointing to theprolate hypothesis.Asasecondtest,Weibelpoints out that, in the case of prolate spheroids, we expect thenearcircularprofilestobemorenumerousthanthe highlyellipticalprofiles.Thiscanbecheckedusingthe axis ratio, r, of an elliptical profile, which is defined as the ratio of the largest diameter to the smallest. In our case, for the sectioned data, the median value of ris 1.52. When this is compared with the minimum (1.02) and maximum (5.70) values of r, the test again confirmstheprolatehypothesis. As discussed in the introduction, the motivation for this research derives from the desire to analyze conformationalchanges,particularly changesin shape ofmitochondriaundergoingapoptosis.Themotivation is to construct a tool for quantifying the progression of shapes after the tubular mitochondrion splits into fragments which eventually become spherical. Our method could be used to monitor those changes quantitatively. We apply our method to a group of HeLacellsin H 2O2-inducedapoptosis. As it turns out, the “proof of concept” presented here is not an idle exercise. Some remarks are in order concerning the a priorifeasibility of the task that we have set out. There is a sense in which the mathematical problem we have formulated has a fundamental ambiguity that may rear its ugly head in certain contexts. Consider two extreme populations of prolate spheroids: ( i) a distribution of nearly perfect spheres with their major axes oriented horizontally, and (ii) a distibution of highly elongated spheroids with all of their major axes nearly vertical. In both cases (since our sectioning plane is defined to be horizontal) the 2D sections are nearly perfect circles. There is thus no information in the shape of the sections that can be used to infer the 3D shape of the spheroids.Thisexampleillustratesthetypeofproblem that can arise when trying to unfold 3D information on both shape and orientation from 2D sections. The presence of such an ambiguity is revealed in the behavior of the likelihood function that is computed from the data. The likelihood function could have several peaks that are not well separated and/or that are of roughly the same height, or it may have an ill defined peak. In these cases this method as it is formulated cannot be applied without bringing in additional information about the population of spheroids that would permit a resolution of the ambiguity.For the data that we analyzed, a simple constraint ofα≤5 for the parameter in Φ0(θ|α)was found sufficient to guarantee a single well-defined interior peak for the likelihood functions L1andL2. This constraint rules out any spheroid population with a standard deviation of θless than 25.0 degrees. That is, if more than about 68% of the spheroids are so flattenedthattheiraxesrisenomorethan25.0degrees above horizontal, then the results we have obtained may not be valid. While we consider this highly unlikely, it has to be recognized as a limitation of this method. We generateda setof micrographs of the group of cells described above and in the Ph.D. thesis of one of the authors (Sun, 2007). The outer membranes of a total of 231 mitochondrial sections were traced in these micrographs using the program ImageJ, which ellipticized each trace yielding minor and major axes, from which a semi-minor axis band an eccentricity e could be determined4. Fig. 2 shows a typical electron micrographusedin ourtracings. 1000 nm Fig. 2.Typical electron micrograph of a HeLa cell havingundergoneH 202inducedapoptosis. Sunet al.(2007) contains a description of how the thin sections were prepared. The thickness of the sections, 80 nm, represents about 15% of a typical particle diameter. Any thoroughgoing effort to use our method would need to address the standard stereological problems of “lost caps” and overprojection associated with the finite thickness of the sections (Weibel, 1979; 1980). We did not do so in this study for two reasons. First, we anticipate that 4ImageJ, (Image Processing and Analysis in Java, http://rsb .info.nih.gov/ij/) uses the “best fit ellipse” algorithm to match the moment of inertia of the regiontracedtodetermine anequivalent el lipse. 18 ImageAnalStereol2009;28:11-22 in our case the corrections would have a minor effect on the results. Secondly, the scope of our enterprise did not extend beyond laying out the theoretical apparatus and testing the feasibility of the full-scale unfolding problem of size, shape and orientation of spheroids. Handling the corrections above for nonsphericalparticlesisnotwelldevelopedandwould havebeen,forusinthiseffort,atheoreticaldistraction. We now recall in outline the procedure described inthepreviousthreesectionsforhandlingthedata,and wesummarizetheresultsasappliedtothedataderived from our micrographs. The predicted size distribution Φ1(r)of the prolate spheroids is dependent on the parameters µr,σr, which describe the most likely normal distribution of the semi-minor axis, r. To produce a likelihood in terms of µr,σr, Eqs. 15 and 17 arecombinedto give L1(D1|µr,σr) =231 ∏ j=1/integraldisplay∞ bjbjΦ1(r|µr,σr)/radicalBig r2−b2 jdr(25) in terms of the 231 semi-minor elliptical axes b1,...,b231obtainedfrom themicrographs. Fig. 3 shows the contour of L1that surrounds the 90% confidence region for µr,σr, as well as the point marking the maximum-likelihood values. Fig. 4 displaysthemaximum-likelihooddistributioninblack. The 90% confidence region is suggested by the green region surrounding the black curve. It is the result of superimposingallofthedistributionswhoseparameter values lie in the 90% confidence region shown in Fig. 3. µr [nm]σr [nm] 230235240245250255260404550556065 Fig. 3.The maximum likelihood point and the 90% confidenceregionforthemeanandstandarddeviation ofthe semi-minoraxis,r,ofthe spheroid.0 100 200 300 400 50000.0020.0040.0060.0080.01 r [nm]Φ1(r) Fig. 4. The maximum likelihood spheroid size distribution (black) and the 90% confidence region (green). The predicted spheroid shape and orientation are found in a similar manner. Eqs. 16, 17, 19, and 21 arecombinedtoproducethelikelihood L2(D2|p,q,α), which is then optimized over p,qandα. The parameters pandqdescribe a Beta distribution, whileαdetermines the orientation distribution. By applying Eq. 20 the optimal p,qyield the eccentricity statistics µe0=0.875 andσe0=0.091. Fig. 5 shows the projection of the 90% confidence region. Fig. 6 illustrates the optimal Beta distribution (in black) as well as the continuum of Beta distributions found within the 90%confidenceregion(in green). µe 0σe 0 0.75 0.8 0.85 0.90.040.060.080.10.120.140.160.180.20.22 Fig. 5.The maximum likelihood point and the 90% confidenceregionforthemeanandstandarddeviation ofthespheroideccentricity,e 0. Fig. 6.The maximum likelihood spheroid eccentricity distribution (black) and the 90% confidence region (green). 19 FLYNNDCET AL:Stereologicalunfoldingmethod 0 50 100 15000.0050.010.015 θΦ0(θ) Fig. 7.The maximum likelihood spheroid orientation distribution (black), a uniform distribution (dashed- blue)andthe 90%confidenceregion(green). The predicted orientation distribution Φ0(θ|α) occurred for α= 2.26. This value of αcorrespondsto σθ=32.82degreesandanisotropyindex,fromEq.22, ofIiso=.838.Thepredicteddistributionoforientation is seen in Fig. 7. Also shown for comparison is the distribution corresponding to perfect isotropy. The 90% confidence analysis shows that our measured anisotropyisstatistically significant. The90%confidenceregionsshowninFigs.3and5 weredeterminedbyperformingaparametricbootstrap as described in the previous section. Explicitly, in the case of spheroid size (semi-minor axis, r, Fig. 3) the maximum-likelihood values µ∗=244 nm and σ∗=52.4 nm were taken to determine the model distribution Φ1(r)from Eq. 18 and the associated predicted distribution for elliptical cross-section size , ˆΨ1(b)after normalizing Eq. 15. A Monte Carlo schemetosimulatetheslicingexperimentwasiterated 5000 times ( N= 5000). Each iteration involved sampling 231 ( m= 231) values of bfrom the distribution ˆΨ1(b)and interpreting that sample as data for a new prediction of µandσ. These 5000 predictions were then processed as explained in the previous section. Fig. 8 shows the raw collection of predictionsinthe (µ,σ)-plane,clusteringnicelyabout theoriginalmaximum-likelihoodprediction(shownin the center). The blue curve in Fig. 9 is an empirical plot, for all 5000 iterations, of LLDversus rank in the form of percentile.AsexplainedinSec.6,thevalueof LLD(= 3.3) determines the boundary of the 90% confidence region.Thisis thecontourshownin Fig.3.220 230 240 250 260 270 280303540455055606570 µr [nm]σr [nm] Fig. 8.Graphical depiction of uncertainty in the method, showing predictions of µrandσrobtained from5000simulationsofthesamplingexperiment(231 sections each) based on a population of spheroids with parameters corresponding to the data-based, maximum-likelihoodprediction(shownin center). 0 0.2 0.4 0.6 0.8 1012345678Log−likelihood Difference Confidence Level Fig. 9.The log likelihood difference as a function of the confidencelevel.The dashedredcurveis obtained from the likelihood ratio test. The solid blue curve is obtainedfromthe parametricbootstrap. The second curve in Fig. 9 (red, dashed) is shown for comparison. If we had chosen not to perform the bootstrap, i.e.not to collect statistics on the performance of our method, we could have used the so-called Likelihood Ratio Test to establish confidenceregions(KendallandStuart, 1983;Hilborn andMangel, 1997). This is an asymptotic result(large m) stating that under certain conditions LLD/2 is distributed like a chi-square statistic of one degree of freedom. The associated percentile curve is shown in red. The discrepancy between these two curves illustrates the necessity, in the case of our method, to base uncertainty results on actual performance statistics,ratherthanassumingtheconditionsrequired forthe LikelihoodRatio Test. A similar parametric bootstrap was performed for results on spheroid eccentricity and orientation. In this case the parameter space ( µe0,σe0,α) is three- dimensional and so is the 90% confidence region. 20 ImageAnalStereol2009;28:11-22 Again, 5000 simulations based on the maximum- likelihood parameters were carried out and the LLD values ranked in an empirical curve. This curve fixed the LLDvalue, and, hence, the contour for the boundinglevelsurfaceoftheconfidenceregion.Fig.5 showsthe projectionofthis regionontothe ( µe0,σe0)- plane,thisprojectiongivingthe90%confidenceregion for thoseparameters. The results from our method may be summarized as follows. At the apoptotic state measured, the mitochondria exhibited a mean eccentricity of 0.875 with a standard deviation of 0.091. In addition the semi-minor axis, r, of the spheroid had a mean value of 244 nm with a standard deviation of 52.4 nm. This corresponds to a semi-major axis of 504 nm. Our most likely prolate spheroid is depicted in Fig. 10. Finally, theprincipalaxesofthe mitochondriashowed a statistically significant directional anisotropy. Their directions showedan isotropicindexof0.838.Thatis, their spreadaboutthe horizontalwas only84% ofthat expectedin anisotropic distribution. Fig. 10.Spheroid with mean dimensions predicted by the method: semi-minor axis 244 nm, semi-major axis 504nm, eccentricity0.875. CONCLUSION We developed a novel stereological method that utilizes a maximum-likelihood scheme to infer the distribution of spheroidal shapes, sizes and orientationsfromtheobserveddistributionofelliptical sections. The maximum-likelihood method is applied to a parametrized set of distributions appropriate for studying the shape of mitochondria undergoing apoptosis. A parametric bootstrap for estimating confidence intervals for the parameters in our distributions is described and illustrated for an example. We expect that future applications to resolving time variations of apoptotic mitochondria will provide insights into apoptosis. We also expect our method, with possible modifications of the parametrization,toserveasausefultemplateforother specializedapplicationsofthe unfoldingproblem.ACKNOWLEDGMENTS This work is the result of interdisciplinary collaboration within the biological, mathematical, computational, and physics departments at San Diego State University. We give special thanks to all the members of the “Mito Research Group”, whose support made this study possible. In particular, we wouldliketothank,JoeMahaffy,PeterBlomgren,and BarbaraA. Bailey. One of the authors (DCF) was supported as a trainee by NIH Roadmap Initiative (Award Number T90DK07015).Microscopystudiesweresupportedby a Blasker Science and Technology Grant awarded by the SanDiegoFoundationtoTGF. REFERENCES Batschelet E (1981). Circular statistics in biology. Londo n: AcademicPress. Bl¨ odner R, M¨ uhlig P, Nagel W (1984). The comparison by simulation of solutions of Wicksell’s problem. J Microsc135:61-74. Cruz-OriveLM(1976).Particlesize-shapedistributions: the generalspheroidproblem.JMicrosc107:233-53. De Hoff RT (1962). The determination of the size distribution of ellipsiodal particles from measurements made on random plane sections. Trans Metall Soc AIME224:474-7. Efron B, Tibshirani RT (1993). An introduction to the bootstrap.New York:ChapmanandHall. Frank S, Gaume B, Bergman-Leitner ES, Leitner WW, Robert EG, Catez F, Smith CL, Youle RJ (2001). The roles of dynamin-related protein 1, a mediator of mitochondrialfission inapoptosis.DevCell,1:515-25. Franklin JN (1977). Some stereological principles in morphometriccytology.SIAMJApplMath33:267-78. Hilborn R, Mangel M (1997). The Ecological Detective: Confronting Models with Data. Princeton, NJ: PrincetonUniversityPress. Keiding N and Jensen ST (1972). Maximum likelihood estimation of the size distribution of liver cell nuclei from the observed distribution in a plane section. Biometrics28:813-829. Kendall M and Stuart A (1983). The Advanced Theory of Statistics, Vol. 2.London:C. Griffin. Mardia KV (1972). Statistics of directional data. London: AcademicPress. Sato E, Kondo N, Wakai F (1996). Particle size, shape and ori- entation distributions: general spheroid problem and application to deformed Si 3N4microstructures. PhilosMagA 74:215-28. 21 FLYNNDCET AL:Stereologicalunfoldingmethod SelbyB(1964).Girdledistributionsonasphere.Biometrik a 51:381-92. Stoyan D, Kendall WS, Mecke J (1987). Stochastic geometry and its applications. Berlin: John Wiley & Sons. SunMG(2007).MitochondrialStructureDuringApoptosis. Ph.D.thesis,San DiegoState University,SanDiego. Sun MG, Williams J, Munoz-Pinedo C, Perkins GA, Brown JM, Ellisman MH, Green DR, Frey TG (2007).Correlated three-dimensionallight and electron microscopy reveals transformation of mitochondria duringapoptosis.NatureCell Biol 9:1057-65.Watson GS (1965). Equatorial distributions on a sphere. Biometrika52:193-201. Weibel ER (1979). Stereological methods, vol. 1. London: AcademicPress. Weibel ER (1980). Stereological methods, vol. 2. London: AcademicPress. Wicksell SD (1925). The corpuscle problem. A mathematicalstudyofabiometricproblem.Biometrika 17:84-99. Wicksell SD (1926). The corpuscle problem. Second Memoir. Case of ellipsoidal corpuscles. Biometrika 18:152-72. 22