Home

Search

Collections

Journals

About

Contact us

My IOPscience

Automated quantification of neurite outgrowth orientation distributions on patterned surfaces

This content has been downloaded from IOPscience. Please scroll down to see the full text. 2014 J. Neural Eng. 11 046006 (http://iopscience.iop.org/1741-2552/11/4/046006) View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 198.91.37.2 This content was downloaded on 22/02/2015 at 02:19

Please note that terms and conditions apply.

Journal of Neural Engineering J. Neural Eng. 11 (2014) 046006 (12pp)

doi:10.1088/1741-2560/11/4/046006

Automated quantification of neurite outgrowth orientation distributions on patterned surfaces Matthew Payne 1,2 , Dadong Wang 1 , Catriona M Sinclair 3 , Robert M I Kapsa 3,4,5 , Anita F Quigley 3,4,5 , Gordon G Wallace 4 , 2,8,9 ¨ Joselito M Razal 4,6 , Ray H Baughman 7 , Gerald Munch and Pascal Vallotton 1 1

Quantitative Imaging Group, Division of Computational Informatics, CSIRO, North Ryde, NSW 1670, Australia 2 School of Medicine, University of Western Sydney, Campbelltown, NSW 2560, Australia 3 Centre for Clinical Neuroscience and Neurological Research, St Vincent’s Hospital, Melbourne, VIC 3065, Australia 4 ARC Centre of Excellance for Electromaterials Science and Intelligent Polymer Research Institute, AIIM Facility, Innovation Campus, University of Wollongong, Squires Way, Fairy Meadow, NSW 2519, Australia 5 Department of Medicine, University of Melbourne, Fitzroy, VIC 3065, Australia 6 Institute for Frontier Materials, Deakin University, Geelong, VIC 3220, Australia 7 Alan G MacDiarmid NanoTech Institute, University of Texas at Dallas, USA 8 Molecular Medicine Research Group, University of Western Sydney, Campbelltown, NSW 2560, Australia 9 Centre for Complementary Medicine Research, University of Western Sydney, Campbelltown, NSW 2560, Australia E-mail: [email protected] Received 11 September 2013, revised 4 March 2014 Accepted for publication 13 March 2014 Published 12 June 2014 Abstract

Objective. We have developed an image analysis methodology for quantifying the anisotropy of neuronal projections on patterned substrates. Approach. Our method is based on the fitting of smoothing splines to the digital traces produced using a non-maximum suppression technique. This enables precise estimates of the local tangents uniformly along the neurite length, and leads to unbiased orientation distributions suitable for objectively assessing the anisotropy induced by tailored surfaces. Main results. In our application, we demonstrate that carbon nanotubes arrayed in parallel bundles over gold surfaces induce a considerable neurite anisotropy; a result which is relevant for regenerative medicine. Significance. Our pipeline is generally applicable to the study of fibrous materials on 2D surfaces and should also find applications in the study of DNA, microtubules, and other polymeric materials. Keywords: neurite outgrowth, neurite anisotropy, orientation distribution, micro-patterned surfaces, image analysis S Online supplementary data available from stacks.iop.org/JNE/11/046006/mmedia (Some figures may appear in colour only in the online journal)

1741-2560/14/046006+12$33.00

1

© 2014 IOP Publishing Ltd

Printed in the UK

M Payne et al

J. Neural Eng. 11 (2014) 046006

of linear structures in two dimensions. We apply it here on a number of image datasets including cells on patterned substrates as well as synthetic images. Our solution was implemented using MATLAB 2010a (The Mathworks, Natick) code, which we have made available as supplementary material (available from stacks.iop.org/JNE/11/046006/mmedia) for this manuscript (a free version of HCA-Vision may also be download at www.hca-vision.com).

1. Introduction Spinal cord injury is the most common cause of paralysis and it affects both the old and the young indiscriminately. The prospects for a potential cure have been considerably improved by the discovery that neuron populations can regenerate under favourable conditions [1]. The main challenge remains to establish the correct connectivity of neurites grown over very long distances. Attention has thus focused on how to guide neurite growth by providing cues—chemical or physical [2–8]. Here, we develop image analysis techniques to characterize the effect of carbon nanotubes on PC12 cells in culture, and we show that the substrate anisotropy provided by the nanotubes translates into a marked anisotropy of their neurite-like extensions. There were significant challenges associated with this exercise, including all the traditional difficulties of neurite tracing itself. The contrast of neurites is typically limited and neurites may feature variable widths and lengths. We have previously addressed such issues and embedded our solution within our software for neurite analysis: HCA-Vision [9, 10]. With regard to the quantification of anisotropy, it is important to include in the analysis only those neurites that are linked to a healthy cell. It is also preferable to retain information on anisotropy at several branching levels, as enabled by HCA-Vision, and it is important that results not be biased by the presence of debris of various nature. The main technical challenge, with regard to the computation of the anisotropy, is that a trace in a digital image grid can present only two orientations (vertical, or horizontal) at every pixel. This information is too coarse-grained to be useful, or even relevant. Thus we fit the digital traces using smooth approximating splines to analytically describe each neurite. Extracting the derivative information from each spline allows us to process the orientation distribution of each neurite on an individual basis. We describe how to choose spline parameters for best results. A survey of the literature shows that existing methods were not suitable for our application. For example, methods involving the Radon [11] and Fourier [8, 12, 13] transforms are affected by variations in neurite contrast and by the presence of cell bodies and debris. Other methods such as estimating the orientation based on the initial or average trajectory [7] are also not satisfactory given that we are dealing with long, flexible neurites. We also aim to do better than regular [5, 2] or random [6] sampling of the orientation along the neurite’s length. To our knowledge, the most applicable and automated tool for measuring the orientation distribution of neurite growth is the ImageJ plugin, OrientationJ [14]. This tensor-based method was shown to be useful for measuring collagen fibre orientation. It builds a histogram of orientations from pixels contained within calculated regions of interest in the image. Unfortunately, results are affected by variations in contrast and thickness and by other objects immediately adjacent to the neurites. Our fully automated methodology overcomes the limitations described above and provides a satisfactory solution to the general problem of quantifying the anisotropy

2. Methods 2.1. Cell culture and staining

PC12 cells were seeded onto carbon nanotube-coated R R plates and plain gold Mylar plates at a gold Mylar −2 concentration of 1200 cells cm , and cultured in Dulbecco’s modified eagle’s medium (DMEM), supplemented with 10% horse serum, 5% fetal bovine serum, 2 mM L-glutamine, 100 U mL−1 penicillin, and 0.1 mg mL−1 streptomycin. The cells were incubated at 5% CO2 and 37 ◦ C for 24 h. The media was then replaced with DMEM supplemented with glutamine/penicillin/streptomycin, with 1% horse serum and 100 ng mL−1 nerve growth factor. After four days of differentiation, the cells were fixed in equal parts of acetone and methanol solution for 10 min at 4 ◦ C, washed in phosphate buffered saline (PBS), and blocked for 1 h at room temperature in PBS with 10% donkey serum. The cells were then incubated with primary antibodies (mouse anti-βIII tubulin and rabbit anti-growth associated protein 43, GAP43) for 2 h and then secondary antibodies (donkey anti-mouse immunoglobulin-594 Alexa Fluor and donkey antirabbit immunoglobulin-488 Alexa Fluor) for 1 h in the same solution, washing with PBS after each incubation. Cell nuclei were then labelled with 4,6-diamidino-2-phenylindole (DAPI), washed in PBS and mounted with Fluoromount (Sigma). A total of 60 images were captured for the carbon nanotube R plates. The antibody plates and 149 for the gold Mylar stains for βIII tubulin and GAP43, together with the DAPI nuclear stain, made up the three colour channels in each image as shown in figures 1(a)–(b). 2.2. Synthetic datasets

We devised four groups of synthetic images representative of different neurite morphologies and anisotropic constraints. Each group contained 50 images of 1600 × 1200 pixels in resolution, similar to those shown in figures 1(c)–(f). The elements pertaining to these images model individual neurite segments, as opposed to complex neuronal structures. As such, these simple lines and curves, or part thereof, can be considered as basic building blocks to produce traces of comparable curvatures and anisotropies. Ground truth orientation distributions used to build the synthetic images were exploited to validate the accuracy of the distributions derived from our image analysis pipeline. The first two image groups depict curved neurites. These groups were constructed using 50 randomly located elliptical (group 1, figure 1(c)) and circular (group 2, figure 1(d)) objects, respectively. The major radii of the ellipses and the radii of 2

M Payne et al

J. Neural Eng. 11 (2014) 046006

(a)

(b)

(c)

(d)

(e)

(f)

R Figure 1. Images of PC12 cells grown on (a) carbon nanotube-coated substrates and (b) plain gold Mylar substrates, along with synthetic

images of (c) group 1: vertically oriented ellipses with an aspect ratio of 0.5 and uniformly distributed major radii, (d) group 2: circles with uniformly distributed radii, (e) group 3: lines with normally distributed orientations and exponentially distributed lengths, and (f) group 4: lines with uniformly distributed orientation and exponentially distributed lengths.

the circles were uniformly distributed between 10 and 100 pixels. The ellipses were oriented in the vertical direction with an aspect ratio of 0.5. Analytical calculations for the ground truth distribution for the elliptical group can be found in the appendix. Clearly, the ground truth distribution for the circle dataset is a uniform distribution. The third and fourth groups contained 100 randomly located linear segments in each image. The orientations of the lines in group 3 were normally distributed about the vertical axis with 10◦ standard deviation, whereas the lines in group 4 were uniformly distributed in all directions. The length

of all lines in both groups was randomly sampled from the exponential distribution: f (x, λ) = λ e−λ(x−xmin ) ,

x  xmin

(1)

with λ = 0.07, to mimic the length distribution displayed by neurite segments on the tailored surfaces (section 2.1), and minimal line length xmin was set to 10 pixels. The thickness of all lines was set to 1 pixel. The ground truth distribution for each group was calculated by summing all length values of similar orientations. To reduce anomalies in the segmentation, the placement of objects in each image prevented overlapping. A single pixel on 3

M Payne et al

J. Neural Eng. 11 (2014) 046006

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 2. A small section from (a) an original colour image of cells grown on CNT—-note the anisotropic growth of the neurite-like xtensions, (b) the enhanced image using steerable filters, (c)–(e) results from HCA-Vision, and (f)–(j) similar images of cells grown ongold R Mylar .

biologically relevant measures and statistics. Further specifics are given in [9] and [10]. From HCA-Vision we obtain skeletonised results of the neurites which are colour coded by (1) the segment identifier, (2) their branching level in the neuronal tree (e.g. primary, secondary, tertiary, etc) and (3) the particular cell they are associated with. This is shown in figures 2(c)–(e). In these results, a neurite segment is defined as the set of connected pixels between branching points. We refer to each neurite segment simply as a neurite in the remainder of this manuscript. To expose the relevant neurite information, we imported the HCA-Vision results into MATLAB. The imported images purposely did not comprise cell bodies and orphan neurites.

each object was also designated as the cell body as HCA-Vision requires this information. The neurite information was in the red channel and the cell body information in the blue channel of the image. All images were smoothed with a Gaussian filter of 5 pixels and were corrupted with random noise of mean zero and 3% variance in all channels. 2.3. Image preprocessing

To aid the tracing of neurites, we first filtered all images using the steerable filter technique described in [15]. Within this setting, we used the three second order derivatives of a two-dimensional Gaussian profile as our steerable filters. By visual inspection, we deemed that a Gaussian width of σ = 1.4 μm achieved the maximum neurite enhancement. The image enhancements are clearly visible in figure 2(b). This step significantly increased the contrast of faint neurites and suppressed background noise. Moreover, images acquired under diverse experimental conditions became consistent in appearance. This permitted analysis of all images (with or without carbon nanotubes) using a single set of image analysis parameters, which is always important if solid biological conclusions are to be drawn. Next, the images were processed by the software HCA-Vision to obtain neurite traces, using its batch processing capability. HCA-Vision is a powerful toolkit designed for high content image analysis and particularly suitable for the study of neuron outgrowth. Its segmentation procedures are based on a non-maximum suppression method [16] and it produces many

2.4. Spline fitting

We encoded the digital neurite traces from each HCA-Vision image as sets of vectors. A single neurite trace consisted of a pair of vectors, x and y, listing the sequence of ordered integer values for the successive pixel coordinates. The MATLAB script ‘splinefit.m’ [17] was applied to these points to produce smooth continuous splines, denoted X(t) and Y(t), where each spline is a set of piecewise cubic polynomials. The splinefit function requires the input of four parameters: the input coordinate vectors x and y, a vector t to parameterize the position of each input coordinate along the neurite length, and a vector k containing a non-decreasing sequence of knot positions to delineate the domain of each 4

M Payne et al

J. Neural Eng. 11 (2014) 046006

polynomial segment. Initially, the vector t was approximated by the chordal distances between the input points, and the vector k contained only the two parametric positions of the neurite’s extremities. In sections 2.4.1 and 2.4.2 we provide details of how we reparameterized t in terms of arc length and increased the number of knots in k to tune the spline’s degrees of freedom. For the spline fitted to a specific neurite, we denote xi, yi and ti as the ith value listed in x, y and t and the point (X(ti), Y(ti)) as the spline approximation of the input point (xi, yi). We also denote the number of pixels in the neurite as m, the total arc length of the spline as l, the number of knots listed in k as nknots, and the number of intervals between the knots as M = nknots − 1. Initially t = [0, . . . , m − 1] and nknots = 2 to give M = 1. The splinefit tool is based on the fitting of cubic B-splines to the input data. B-splines are built as a linear summation of a set of basis functions: i.e. an nth order B-spline is given by: s(t ) =

M 

di Bni (t )

when solved, determines the least squares solutions for the weighting coefficients of equation (2). When selecting the parameterization scheme, we considered the bias caused by the non-uniformity in the spacing of the input points along the neurite. For example, the distance between √ the estimations of the input points decreases by a factor of 2 as the neurite is rotated by 45◦ from the horizontal or vertical directions (see insert of figure 3(c)). As a result, schemes such as the uniform parameterization would fit the trace more closely in sections that are either vertical or horizontal. We therefore require an arc length parameterization to account for the geometric spacing of the input points. Arc length parameterization has also been shown to be the best parameterization method for parametric interpolation [20], which is utilized extensively in the production of our orientation distributions. It also ensures that knots placed at regular parametric intervals are represented by arcs of equal length, thus distributing the smoothing properties of the spline evenly along the neurite’s length. Arc length parameterization requires the elements of t to range on the interval [0, l] whereby each element ti ∈ t stores the distance traversed along the spline from the origin (X (0), Y (0)) to the point (X (ti ), Y (ti )). Given however that X (t ) and Y (t ) are cubic polynomials, the task is to determine the solution to the arc length integral, as shown below in equation (6).      t  d(X (t )) 2 d(Y (t )) 2 + dt. (6) l= dt dt 0

(2)

i=1

where di is a weighting coefficient determined by input points and Bni (t ) is the ith basis function of order n. The basis functions are only dependent on the sequence of knot values spanning the domain of t. Given the knot sequence k = [k1 , . . . , knknots ], the first order basis functions are defined as:  1, ki  t < ki+1 0 . (3) Bi (t ) = 0, otherwise

To achieve arc length parameterization, we applied a spline refitting process that used a reparameterization of the vector t upon each repetition. The initial parameterization was estimated using a chordal parameterization scheme whereby the elements in t were defined as:

Higher order basis functions are evaluated as a combination of two lower order basis functions: t − ki ki+n − ti n−1 Bni (t ) = Bn−1 (t ) + B (t ). (4) ki+n−1 − ki i ki+n − ki+1 i+1 Consequently, a basis function of order n has support over n + 1 knot intervals and the knot interval [ki, ki + 1] is represented by a polynomial requiring the summation of n + 1 non-zero basis functions. The weighting coefficients for the summation in equation (2) are obtained using the least square method. For example, the following equation is solved to fit X(t) in the least squares sense [18]: ⎤⎛ ⎞ ⎛ ⎞ ⎡ n d1 x1 B1 (t1 ) · · · BnM (t1 ) ⎥ ⎜ .. ⎟ ⎜ .. ⎟ ⎢ .. .. .. (5) ⎦⎝ . ⎠ = ⎝ . ⎠. ⎣ . . . Bn1 (tm )

···

BnM (tm )

dM

ti+1

t1 = 0 = ti + (xi+1 , yi+1 ), (xi , yi )

(7)

where  ·  represents the Euclidean distance between the two points. This initial parameterization formed the splines X0 (t ) and Y0 (t ) using the splinefit tool described previously. We use the reparameterization detailed in the appendix to fit subsequent splines, labelled X j (t ) and Y j (t ) for j = 1, . . . , nrefits . The improvements resulting from the reparameterization are shown in table 1. In table 2 we show the improvements of fitting accuracy and the increase in correlation between the spline parameter and arc length after several iterations of the refitting process applied to the neurite in figure 3. We also see that convergence occurred after the third iteration, thus we used three refitting iterations on our image datasets.

xm

Ultimately, the two key factors controlling the shape of a spline are the parameterization of the vector t and the placement of knots contained within k for each spline. We discuss our optimization methods for these control factors in the following sections. The vector t is a list of nondecreasing values selected to parameterize the locations of the ordered input coordinates contained in the vectors x and y. The choice of parameterization scheme is a major factor in determining the shape of the B-splines [19]. In particular, the elements of the parameterization vector t are utilized in the formation of the matrix shown in equation (5) which,

2.4.1. Spline parameterization.

The other major factor affecting the shape of the splines is the positioning of the knots in the vector k. The basis functions used to construct the spline are solely dependent on the spacing of the knots. Also, a cubic polynomial is assigned to each knot interval, thus defining the spline’s degrees of freedom. The ideal spline contains

2.4.2. Knot placement.

5

M Payne et al

J. Neural Eng. 11 (2014) 046006

(a)

(c)

(b)

Figure 3. An example of the spline fitting process to a neurite segment with m = 196 pixels using M = 4 knot intervals. The splines X (t )

and Y (t ) are fitted to the x and y pixel coordinate vectors and are shown as functions of t in (a) and (b) respectively. The splines consist of cubic polynomial segments (alternating between red and green at each segment) with each segment terminating at the knots listed in the vector k (shown in blue). The resulting smoothed neurite trace is shown in (c) and measures a length of li = 146.58 pixels. The arrow indicates the direction of increasing t. The inserts shows the mapping of the input coordinates at√different locations of the neurite to their respective estimation on the smoothed spline. The distance d1 at 0◦ and 90◦ will on average be 2 longer than d2 at 45◦ and 135◦ . Thus arc length parameterization is required for correct sampling. Table 1. An outline of our method described in section 2.4 to adjust the factors controlling the shape of a spline. The loop from lines 7 to 20 calculates the arc length parameterization obtained in section 2.4.1 whilst the loop from lines 3 to 22 increments the number of uniformly spaced knots required in section 2.4.2.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

Set nknots = 1 Set t as the cumulative sum of chordal lengths between {(xi , yi )|i = 1, . . . , m} DO Increment nknots Set k uniformly Fit splines X0 (t ) and Y0 (t ) using x, y, t & k FOR j = 0 TO number of refits - 1 Calculate X j (t ), Y j (t ), X j (t ) & Y j (t ) FOR number of Newton’s method iterations FOR i = 1 TO m xi = X j (ti ) − xi yi = Y j (ti ) − yi fi = xi X j (ti ) + yiY j (ti ) fi = X j (ti )2 + Y j (ti )2 + xi X j (ti ) + yiY j (ti ) Update ti with ti − fi / fi NEXT NEXT Set t as the cumulative sum of chordal lengths between {(X j (ti ), Y j (ti ))|i = 1, . . . , m} Fit splines X j+1 (t ) & Y j+1 (t ) using x, y, t & k NEXT Find maximum residual between the splines and the input points UNTIL maximum residual < tolerance

6

M Payne et al

J. Neural Eng. 11 (2014) 046006

Table 2. The effect of spline reparameterization applied to the neurite segment in figure 3. The average and maximum residuals are the respective distances (in pixels) of the input coordinates to their respective estimation on the spline. These columns indicate how closely the spline is fitting the input data. The average and variance in arc length per parametric unit show the improvements in our arc length parameterization on each iteration. The arc length of the neurite was based on the length of polyline consisting of 1000 nodes per parametric unit.

Iteration

Average residual

Maximum residual

Average arc length per parametric unit

Variance of arc length per parametric unit ( × 10−3)

0 1 2 3 4 5

0.7165 0.5784 0.5646 0.5628 0.5626 0.5626

3.9109 2.6907 2.6195 2.6105 2.6094 2.6093

0.7446 0.9961 0.9941 0.9938 0.9938 0.9938

68.1366 2.7876 2.5587 2.5481 2.5481 2.5482

enough knots to fit closely to the input points (x and y) whilst moderating oscillations in the data by limiting the degrees of freedom. This is vital in obtaining the true orientation of the neurites. To obtain the appropriate knot sequence for a specific neurite we assess the fitting accuracy of splines fitted with increasing numbers of uniformly spaced knots. Starting with a spline consisting of a single knot interval, we retrieve the maximum residual:    rmax = (8) max xi2 + y2i .

parametric locations of points along the spline that produced an orientation belonging to the set:   π 2b − 1 , b = 1, . . . , nbins (9) = 2 nbins where nbins was the number of uniformly spaced divisions in the orientation distribution over the range [0, π ]. We used nbins = 180 to produce a distribution with 1◦ intervals. For a specific neurite segment, we defined the vectors T and θ to record the parametric positions and orientations evaluated at these partition points, i.e.: ⎡ ⎤ ⎡ ⎤ θ1 T1 ⎢θ2 ⎥ ⎢T2 ⎥ (10) T = ⎣ ⎦, θ = ⎣ ⎦, .. .. . .

i=1,...,m

This approximates the maximum distance between the spline and the control points. If the maximum residual exceeds a user specified tolerance, the number of uniformly spaced knots is incremented and the process is repeated until the condition is satisfied. The tolerance needs to consider the image resolution, the complexity of the neurites’ traces, and the signal to noise ratio of the images. A smaller tolerance has a higher computational demand; with a perfectly fitting spline (rmax = 0) achieved when using m − n knots. However, given a reasonable tolerance for adequate smoothing, a solution is typically found relatively quickly. We found that a tolerance of 3 pixels (i.e. 2.24 μm in our particular images) provided a good balance between accuracy and smoothing for splines fitted to our images. Complex neuronal trees exhibiting long neurite segments and frequent sharp directional changes may need to consider non-uniform rational B-splines (NURBS) for more accurate curve fitting. There are several methods to automate the placement of non-uniformly distributed knots in this generalization, for example [21] and [22]. We found however that 97.8% of all neurites in all images satisfied the previous tolerance after the first iteration, such that NURBS were not considered further.

with the relationship:



 X  (Ti ) ∈ (11) Y  (Ti ) where X  (t ) and Y  (t ) are the derivative functions belonging to the spline. We calculated the values of T by processing the spline derivative coefficients in the form of large matrices—the processing was performed under MATLAB, an environment where vectorized calculations are performed much faster than loops. First we created the matrices CX and CY to contain the coefficients from X (t ) and Y (t ) respectively. For example the M cubic polynomials in X (t ) for a particular neurite make up the matrix: ⎤ ⎡ 3a1 2b1 c1 ⎢ .. .. ⎥ (12) CX = ⎣ ... . . ⎦ θi = tan−1

3aM

2bM

cM

where the kth segment of the spline is a polynomial of the form: X(k) (t ) = ak t 3 + bk t 2 + ck t + dk .

2.5. Fast derivative processing

(13)

Vectors containing the upper and lower knot positions for the polynomials described by each row in CX and CY were also recorded. The CX and CY matrices were combined into a single matrix: ⎤ ⎡ C1 ⎥ ⎢ (14) C = ⎣ ... ⎦

Describing the neurite segments as analytic cubic functions permits the precise orientation of the neurite to be found at any position along its entire length. Moreover, to apportion a spline’s length to a discrete number of intervals (i.e. histogram bins in the orientation distribution) the splines were partitioned at points yielding an orientation equal to any of the distribution’s interval boundaries. That is, we identified the

Cnbins 7

M Payne et al

J. Neural Eng. 11 (2014) 046006

where each submatrix Cb utilized a different θb ∈  in the form: Cb = (tan(θb )CX − CY ).

pixel values. That is, all lines within a small angular range are represented by a single pixel configuration. The resolution angle is dependent on length and is greatest at 0◦ and 90◦ . Furthermore, it can be significantly larger than the bin width of the orientation distribution causing a portion of the spline’s length that spans several bins to be clumped into a single bin. To remedy the issue, we smoothed each neurite’s orientation distribution vector with a unique smoothing rectangular window element based on its maximum angular resolution error. The width of the smoothing element for a given neurite was calculated as:   2nBins 1 w= (17) tan−1 π l−1 where l is the length of the neurite. This has the effect of ensuring that the length of the s(t ) was uniformly distributed over all bins contained within the angular resolution range instead of concentrated in specific bins (such as 0◦ and 90◦ ). The smoothed distributions for all neurites were summed to give an overall distribution per image. The average distributions for the two tailored surface and four synthetic datasets from sections 2.1 and 2.2 respectively were then computed and normalized in terms of total neurite length. In each normalized distribution the total neurite length is equal to that of a perfectly isotropic distribution producing a unit length on all orientations.

(15)

Each row of the matrix C was treated as a rearranged version of equation (11) and solved to obtain the required solutions for T. The solutions were found by applying column operations to C in accordance with the quadratic formula to obtain two solutions per row (rows with a zero element in the first column were treated separately to produce only a single solution). Only real solutions complying with their recorded upper and lower knot boundaries were sorted and stored in the vector T. The corresponding orientation θb for each solution in T was stored at the same index in the vector θ. We also appended the parametric positions and orientations of the neurite extremities (t = 0 and t = l) to their respective vectors. A further processing speedup was obtained by concatenating the matrix C in equation (14) for all splines in the image into a single matrix. This processed all splines in an image in a single iteration. An additional vector was used to record the spline identifier for each row in the composite matrix such that the solutions could be decomposed into their individual T and θ for each neurite. 2.6. Generating the orientation distribution

The orientation distribution was defined to have divisions at the orientations contained in  as defined in equation (9), causing each bin to be centred at the orientations: bπ φb = , b = 1, . . . , nBins . (16) nBins To allocate a neurite’s length to the correct bins, the positions listed in T were used to partition the curve (X (t ), Y (t )) into a number of small arcs. We labelled si as the arc limited by Ti and Ti+1 . The arc length parameterization allowed us to calculate the length of si as |Ti+1 − Ti |. As the splines were continuous and smooth, the orientation along the arc si remained between the orientations θi and θi+1 . Thus the length of si was assigned to the bin interval φb located between θi and θi+1 . The only exception was when θi and θi+1 were equal. In this case, we chose the φb closest to the orientation calculated at the midpoint of the arc. The orientation distribution for each neurite was stored separately. The total area below each distribution is the length of its corresponding neurite. The orientation distribution for each image was simply calculated as the summation of the distributions calculated for each individual neurite contained within the image. We found, however, that this could lead to the formation of peaks in the horizontal (0◦ ) and vertical directions (90◦ ) in cases containing straight neurites close to these orientations. Similar artefacts were also visible in the results obtained by a comparable quantification tool: OrientationJ. Even the results produced from the same set of images rotated by an arbitrary angle presented similar peaks at 0◦ and 90◦ (as opposed to these peaks shifting together with the rotation angle, which would be expected to happen if these peaks were genuine). We attributed the cause of these peaks to an angular resolution limitation when trying to represent short, straight lines using discrete

3. Results Only neurites longer than 7.5 μm (10 pixels) were processed to ensure that only established neurites subjected to guidance cues (when present), were included in our results. Results for the six image datasets are displayed in figure 4. We compared our distributions with those produced by OrientationJ [14]. OrientationJ is a plug-in for the popular image analysis tool ImageJ, developed to measure the local angular distribution of collagen fibres in 2D images. The methodology is based on the calculation of local orientation, energy and coherence from the Hessian matrix at every pixel in the image. Energy and coherence thresholds generate an image mask of significant regions in which a histogram of orientations is built. For our synthetic images we applied energy and coherence thresholds equal to 30% and 70% respectively. For the patterned substrate images, these thresholds were reduced to 1% and 50% to include faint neurites in the analysis. In all cases, the image derivatives were evaluated using cubic spline gradient interpolation. For comparison purposes, the OrientationJ results are also normalized in a similar manner to section 2.6. The full complement of results from our method and OrientationJ for all groups can be seen in figure 4. Our method produced a significantly larger peak at 90◦ orientation in figure 4(a). To objectively evaluate the accuracy of both methods, ground truths were established for all synthetic image groups. The ground truth distribution for all ellipses in group 1 were calculated by numerically evaluating the elliptic integral as done in the appendix. We assumed a uniform distribution for the circles in group 2. For the anisotropic and isotropic lines in 8

M Payne et al

J. Neural Eng. 11 (2014) 046006

Figure 4. Normalized results from our spline-based method and OrientationJ (ground truths also shown where available) for (a) 60 images R substrates, and 50 synthetically of PC12 cells grown on carbon nanotube substrates, (b) 149 images of PC12 cells grown on gold Mylar constructed images consisting of (c) ellipses (aspect ratio of 0.5), (d) circles, (e) lines orientated with normally distributed orientations and (f) lines orientated with uniformly distributed orientations.

our method in general, can more accurately find the orientation of the neurite centrelines. Applying circular statistics [7] to the images from patterned substrates, we obtained a mean of 91.6◦ and standard deviation of 21.9◦ for neurites grown on carbon nanotubes. For R substrates we obtained a standard the isotropic gold Mylar ◦ deviation of 76.3 .

Table 3. The total sum of squares with respect to the ground truth for each group of synthetic images.

OrientationJ Our method

Ellipses

Circles

Anisotropic lines

Isotropic lines

26.1 2.0

27.7 0.008

72.3 20.8

13.4 7.7

4. Discussion

group 3 and group 4, we calculated weighted histograms using the length and orientation data stored for each line during the construction of the images. The ground truth distributions are also shown in figure 4. To quantify the accuracy of the distributions produced by our method and OrientationJ, we calculated the sum of squares for the difference between each method and the ground truth. This data is available in table 3. All four cases resulted in an error less than that achieved by OrientationJ, indicating that

Our measurement of anisotropy on the basis of discrete neurite traces rather than on the basis of averages formed over entire regions of interest (as done e.g. by OrientationJ) affords several advantages. For example, in the extreme case of very faint neurites, one may imagine that the variable anisotropy along the neurites may lead to low anisotropy values, even if all neurites were in reality oriented parallel to each other. Using HCA-Vision, we were able to detect the presence of very faint 9

M Payne et al

J. Neural Eng. 11 (2014) 046006

Figure 5. Our results for the PC12 cells grown on CNT (figure 4(a)) with the orientation distribution of each neurite-like segment ordered by

length before being summed and normalized. We see that longer segments display greater anisotropy in the direction of the substrate’s guidance cues.

element related to its length. Although alternative means exist to correct for these anomalies [23], our simple correction technique significantly reduced the presence of artificial peaks within our distributions. Prior versions of our method suffered from the same symptoms until we implemented the techniques described in sections 2.4.1, 2.4.2 and 2.6. This issue was such a challenge that in an early version of our algorithm, we systematically rotated the image and bundled the corresponding results in an effort to limit these effects. To eliminate this step, we experimented with different knot intervals and parameterization methods to improve the efficacy of the spline fitting process. This significantly reduced the presence of the abnormalities, particularly as the knot interval increased. The remaining artefacts were further reduced as we implemented the smoothing process from section 2.6. This smoothing was designed to have a stronger influence on short straight neurites whereby the spline fitting process has limited flexibility in its parameters. Thus, our pipeline presented in this contribution avoids the image rotation step and is, therefore, considerably more efficient. Our method adds only few parameters to those already necessary to generate traces under HCA-Vision. The main parameter represents the maximum distance from the spline curve to pixels belonging to the discrete traces. It is difficult to specify precise guidelines for choosing this parameter as it is influenced by our own expectation of neurite smoothness formed by examining a particular data set. We recommend starting with the default value of 3 pixels and to modify it only if necessary. For example, a controlled surface that would induce sharp bends in neurites would require lowering the tolerance—otherwise, the bends would possibly be smoothed out by the spline. We hope to experiment with such complex surfaces in the future.

neurites and represent their centrelines in an unbiased fashion. This is necessary to avoid over-representing the brighter and thicker neurites in the orientation distributions. HCA-Vision only considers neurites connected to cell bodies—an important aspect in removing the influence of debris, of cell bodies, and of miscellaneous noise contained in the image. An individual treatment of neurites additionally provided the means to correlate orientation distribution with other neurite measures. For example, figure 5 illustrates the decomposition of our result shown in figure 4(a) into distributions based on individual neurite length. This demonstrates that longer neurites are more susceptible to the guidance cues generated by the CNT substrate. The process of fitting splines to the neurite traces tends to remove any unwanted oscillations due to noise. These analytic representations of the neurites allowed the precise orientation to be calculated at any location along the neurite, which was essential for the accurate partitioning of the spline’s length into finely graded intervals in the orientation distribution. Furthermore, the arc length parameterization gave an efficient means to calculate the length of each section and to ensure that all orientations were represented in an unbiased fashion. This process is of particular benefit when dealing with curved neurites as seen by the lower errors attained in group 1 and group 2 in table 3. Our results on synthetic data sets have also revealed the presence of spurious oscillations at periodic intervals when using OrientationJ instead of our methodology. It was difficult to trace the origin of these oscillations, but it is probable that they stem from the resolution limitations imposed by the discreteness of the image raster. We have taken particular care to minimize these effects in our results through arc length parameterization, knot placement and by smoothing each neurite’s orientation distributions with an 10

M Payne et al

J. Neural Eng. 11 (2014) 046006

The integral of equation (A.5) is in the form of an elliptic integral of the first kind, thus we apply a simple application of the trapezoidal rule to find a numerical approximation. The solution for π2  θ < 2π is obtained by symmetry.

5. Conclusion In the field of cell-based assay, there has been a sustained effort to move from average measurements performed over entire images, to ‘per cell’ measurements. Only then may the cell-to-cell biological variability be adequately appreciated. Accordingly, biologists typically produce cell-based statistics when they evaluate their images manually. Our neurite anisotropy measurement method is in line with these requirements. The present work was made possible by the availability of neurite analysis software, such as HCA-Vision, which delivers neurite information on a per cell basis. In the near future, we will endeavour to integrate the anisotropy functions described here into HCA-Vision. This will allow for the formulation of database queries involving anisotropy in addition to standard features such as neurite length, cell body diameter, neurite thickness, etc. For example, it will be possible to investigate whether thinner neurites tend to orient better along guidance cues. In the field of drug screening, there is considerable interest in compounds that modulate the anisotropy of the cytoskeleton, whether actin filament, microtubules, or intermediate filaments, and our new tool should help there too.

Arc length reparameterization

At each iteration we calculated the first and second derivatives of the spline, namely X j (t ), Y j (t ), X j (t ) and, Y j (t ), and computed the vectors x and y, defined as: x = X j (t) − x y = Y j (t) − y.

The least squares solutions of equation (5) define the splines X j (t ) and Y j (t ) such that the sum of the squares of all elements in the x and y vectors were minimized. This could be considered the best global fit of the input data in x and y, given the current parameterization in t. However, given that t is only an approximation which we are trying to refine, there may be better choices of t that will minimize each element in x and y further than that obtained in equation (A.6). Thus, we refined the estimation of t by adjusting each element ti ∈ t such that the coordinate (X j (ti ), Y j (ti )) on the spline was at a minimal distance from the input point (xi , yi ). That is, for each element pair xi ∈ x and yi ∈ y, we solved:  d 2 xi + y2i = 0. (A.7) fi = dt The solution was achieved by applying Newton’s method in the form: fi ([ti ]k ) (A.8) [ti ]k+1 = [ti ]k −  fi ([ti ]k ) where: fi (t ) = xi X j (t ) + yiY j (t )

Appendix Estimation of ellipse orientation distribution

An ellipse was described by the following equations: dx = − sin θ dθ dy y = α sin θ = α cos θ (A.1) dθ where α is the aspect ratio. The arc length element on the ellipse was:  ds = dx2 + dy2  = sin2 θ + α 2 cos2 θ dθ (A.2)  2 2 = 1 + (α − 1) cos θ dθ x = cos θ

and its orientation (for 0  θ < π2 ) as:   −dx π + φ = tan−1 dy 2   sin θ π + . = tan−1 α cos θ 2 Rearranging and deriving gives: θ = tan−1 (−α cot(φ)) α dθ = . 2 dφ 1 + (α − 1) cos2 φ Thus, ds ds dθ = · dφ dθ dφ =

1+

(α 2

α − 1) cos2 φ

.Z.

(A.6)

fi (t ) = X j (t )2 + Y j (t )2 + xi X j (t ) + yiY j (t ) (A.9) and [ti ]k represents values of ti after the kth iteration of Newton’s method. In practice, this method provided fast convergence and, as only relatively small adjustments were required from our initial estimation, convergence issues did not occur. For our images we conservatively fixed the number of Newton iterations to 5. To achieve the arc length parameterization given in equation (6) the values of t must obey: t =0  ti+1 1 (A.10) ti+1 = ti + X j (t )2 + Y j (t )2 dt. ti

However, given our efforts to accurately refine the positions of each element in t, we approximated the integral in equation (A.10) using the chordal distances between the estimations found on the spline, i.e.: t¯1 = 0 t¯i+1 = t¯i + ||(X j (ti+1 ), Y j (ti )), (X j (ti ), Y j (ti ))|| (A.11) where t¯i is the arc length reparameterization of ti . The high frequency of points on the spline was sufficient to portray a polyline closely matching this spline, and thus close to the optimal parameterization [20]. The new parameterization in t¯ was used to refit the next iteration of splines X j+1 (t ) and Y j+1 (t ). The process from equation (A.6) was then repeated for subsequent iterations.

(A.3)

(A.4)

(A.5) 11

M Payne et al

J. Neural Eng. 11 (2014) 046006

[12] Tunak M and Linka A 2007 Analysis of planar anisotropy of fibre systems by using 2D Fourier transform Fibres Text. East Eur. 15 86–90 www.fibtex.lodz.pl/article249.html [13] Redon C, Chermant L, Chermant J L and Coster M 1999 Automatic image analysis and morphology of fibre reinforced concrete Cement Concr. Compos. 21 403–12 [14] Rezakhaniha R, Agianniotis A, Schrauwen J T C, Griffa A, Sage D, Bouten C V C, van de Vosse F N, Unser M and Stergiopulos N 2012 Experimental investigation of collagen waviness and orientation in the arterial adventitia using confocal laser scanning microscopy Biomech. Model. Mechanobiol. 11 461–73 [15] Jacob M and Unser M 2004 Design of steerable filters for feature detection using canny-like criteria IEEE Trans. Pattern Anal. Mach. Intell. 26 1007–19 [16] Sun C and Vallotton P 2009 Fast linear feature detection using multiple directional non-maximum suppression J. Microsc. 234 147–57 [17] Lundgren J 2007 Splinefit Matlab Central, Mathworks www.mathworks.com.au/matlabcentral/fileexchange/ 13812-splinefit [18] Bishop C M 2006 Pattern Recognition and Machine Learning (Information Science and Statistics) (New York: Springer) [19] Haron H, Rehman A, Adi D I S, Lim S P and Saba T 2012 Parameterization method on B-spline curve Math. Probl. Eng. 2012 640472 [20] Kuznetsov E B and Yakimovich A Y 2006 The best parameterization for parametric interpolation J. Comput. Appl. Math. 191 239–45 [21] Jacobson T J and Murphy M J 2011 Optimized knot placement for B-splines in deformable image registration Med. Phys. 38 4579–82 [22] Park H and Lee J H 2007 B-spline curve fitting based on adaptive curve refinement using dominant points Comput. Aided Des. 39 439–51 [23] Nguyen C V, Nguyen T D, Wells J C and Nakayama A 2010 Interfacial PIV to resolve flows in the vicinity of curved surfaces Exp. Fluids 48 577–87

References [1] Reynolds B A and Weiss S 1992 Generation of neurons and astrocytes from isolated cells of the adult mammalian central-nervous-system Science 255 1707–10 [2] Dowell-Mesfin N M, Abdul-Karim M A, Turner A M, Schanz S, Craighead H G, Roysam B, Turner J N and Shain W 2004 Topographically modified surfaces affect orientation and growth of hippocampal neurons J. Neural Eng. 1 78–90 [3] Fan L, Feng C, Zhao W, Qian L, Wang Y and Li Y 2012 Directional neurite outgrowth on superaligned carbon nanotube yarn patterned substrate Nano Lett. 12 3668–73 [4] Jang M J and Nam Y 2012 Geometric effect of cell adhesive polygonal micropatterns on neuritogenesis and axon guidance J. Neural Eng. 9 046019 [5] Kofron C M, Liu Y T, Lopez-Fagundo C Y, Mitchel J A and Hoffman-Kim D 2010 Neurite outgrowth at the biomimetic interface Ann. Biomed. Eng. 38 2210–25 [6] Lee Y S, Collins G and Arinzeh T L 2011 Neurite extension of primary neurons on electrospun piezoelectric scaffolds Acta Biomater. 7 3877–86 [7] Li G N and Hoffman-Kim D 2008 Evaluation of neurite outgrowth anisotropy using a novel application of circular analysis J. Neurosci. Methods 174 202–14 [8] Jang M J, Namgung S, Hong S and Nam Y 2010 Directional neurite growth using carbon nanotube patterned substrates as a biomimetic cue Nanotechnology 21 235102 [9] Wang D D, Lagerstrom R, Sun C M, Bishof L, Valotton P and Gotte M 2010 HCA-vision: automated neurite outgrowth analysis J. Biomol. Screen. 15 1165–70 [10] Vallotton P, Lagerstrom R, Sun C, Buckley M, Wang D, De Silva M, Tan S S and Gunnersen J M 2007 Automated analysis of neurite branching in cultured cortical neurons using HCA-vision Cytometry A 71 889–95 [11] Krause M, Hausherr J and Krenkel W 2011 Computing the fibre orientation from Radon data using local Radon transform Inverse Probl. Imaging 5 879–91

12