Elucidating the control and development of skin patterning in cuttlefish.

Few animals provide a readout that is as objective of their perceptual state as camouflaging cephalopods. Their skin display system includes an extensive array of pigment cells (chromatophores), each expandable by radial muscles controlled by motor neurons. If one could track the individual expansion states of the chromatophores, one would obtain a quantitative description-and potentially even a neural description by proxy-of the perceptual state of the animal in real time. Here we present the use of computational and analytical methods to achieve this in behaving animals, quantifying the states of tens of thousands of chromatophores at sixty frames per second, at single-cell resolution, and over weeks. We infer a statistical hierarchy of motor control, reveal an underlying low-dimensional structure to pattern dynamics and uncover rules that govern the development of skin patterns. This approach provides an objective description of complex perceptual behaviour, and a powerful means to uncover the organizational principles that underlie the function, dynamics and morphogenesis of neural systems.

Cuttlefish and octopus have an unmatched ability to change their external appearance for camouflage or communication1. When camouflaging, they produce a statistical approximation of their visual environment following rules that remain unknown. Because cephalopod camouflage appeared as a response to predators and because their performance can fool humans as well, the rules of pattern generation that they express may be instructive about texture perception across animals, and reveal biological solutions to a general problem of computational vision and neuroscience2–6.

Since the pioneering work of Florey on cephalopod chromatophores7,8, Hanlon, Messenger and colleagues have revealed the remarkable complexity of this system,9–11. Pigment-carrying chromatophores—the pixels of this 2D texture generation system— expand and contract in direct response to the activity of motoneurons8, which project from the brain12 and make excitatory glutamatergic synaptic connections13 onto sets of muscles arranged radially14. Chromatophores operate in concert with other specialized cells (e.g., leucophores and iridophores) and dermal muscular systems to generate a rich array of coordinated textures, dynamic patterns and behaviours.

The rules of neural control governing this system remain largely elusive, owing greatly to the challenges of tracking large numbers (thousands to millions) of small (15-100 µm diameter) chromatophores in soft-bodied, behaving animals. Because each chromatophore is controlled by a small number of motoneurons and conversely, because each motoneuron controls a small number of chromatophores (its motor unit, or MU)12,14,15, we reasoned that chromatophore expansion could serve as a proxy for motoneuron activity. Analysis of the joint statistics of chromatophore variation might in turn reveal the structure of a hypothetical control hierarchy. This study presents our solutions to this challenge, a method for tracking nearly all chromatophores of a cuttlefish’s dorsal mantle at high frame rate and over developmental timescales. Using this technique, we take the first quantitative steps towards elucidating the control, dynamics and morphogenesis of this system.

Tracking chromatophores in free behavior

Freely behaving animals were filmed in a tank with variable backgrounds (Methods). Our first goal was to segment all visible chromatophores in all images and then align these images by mapping them into a common reference frame (Fig. 1a). In any recording session, continuous image sequences with the animal in view and in focus (“chunks”, Fig 1b) were interrupted by unusable periods of movement (grey, Fig. 1b). Chunks were selected post-hoc using a statistic of focus (Methods, Fig. 1b). Within a chunk, the pixels of each frame were first classified16 as belonging either to a chromatophore (of any colour) or to background (Methods, Fig. 1c, Extended Data Fig. 1). All the frames in one chunk were mapped into a common reference frame using sparse optical flow17 (Fig. 1d, Supplementary Video 1) and averaged into one “master frame” (Fig. 1e, top row, left and middle).

Figure 1

Chromatophore tracking in behaving cuttlefish.

a, Schematic of approach to track the expansion state of many single, identified chromatophores (e.g. *) over time through nonlinear image alignment. b, Data acquisition. ‘Chunks’ of in-focus image segments are identified with a focus statistic (black, bottom), separated by out-of-focus segments (grey). c, Segmentation. Pixels are classified as belonging either to a chromatophore (black, right) or to background (white). d, Alignment. Compare the averages of two frames with and without within-chunk alignment. e, Stitching. Top row: one master frame (middle) is mapped into the reference frame of another (left), resulting in image at right. Middle two rows: Corresponding patches of skin (purple and green squares at top), after alignment (left and middle). Matching made possible by peak in spatial cross-correlation function (right) even when chromatophore states differ (green). Bottom row, left: shifting or rotating corresponding patches (purple frame) rapidly decreases correlation (Pearson’s). Bottom right: distribution of correlation values between 64x64 pixel skin patches (as above) across aligned master frames at all translations (sampled every pixel), and rotations (sampled every 60 degrees). Correct matches (red): N=36. Non-matching pairings (black), N=485,722,908. f, Defining chromatophores. Left: Queen frame (N=167,065 aligned, averaged, segmented frames). Right: zoom-in of the red square at left; single-chromatophore centers indicated with white dots. g, Raw expansion state of chromatophore in f (blue) over successive frames.

Extended Data Figure 1

Accuracy of chromatophore classifier

a, Test patch of skin used for classifier testing. b, Segmentation by expert human. c, Segmentation by classification algorithm. d, Composite image comparing manual (annotation) and automatic (prediction) segmentation. 87%-pixel agreement, with differences mostly on the edges surrounding chromatophores. e, Quantification of region overlap. Regions defined from watershedding the composite image shown in d. Correct detection: regions labelled by both methods. False detections: regions identified by automatic but not by manual segmentation. Failed detections: regions identified by manual but not automatic segmentation. f, For all regions, annotated vs. predicted size. Line: identity.

To track individual chromatophores across filming gaps, we stitched chunks together. By correlating a small patch of skin (purple and green frames, Fig. 1e) within one master frame with all possible positions and orientations of identically-sized patches in another master frame, a single “fingerprint” match was usually detected (see 2D correlations in Fig. 1e; note sensitivity to small shifts: bottom left, Extended Data Fig. 2). We used a matching procedure over a grid of patches (Fig. 1e, top left and bottom right, red), interpolating between matching points, to map all master frames into a common reference frame (Supplementary Video 2, Methods). Over 970 master frames analysed, 85% had an average mapping error of ≤ 3 pixels /20±6 µm. The following data come from 1,178,146 frames so mapped, from 6 animals.

Extended Data Figure 2

Sensitivity of correlation between skin patches to small image translations and rotations

Left: Skin patches from the 2 sets of matching master frames shown in Fig. 1e. Middle, right: Composite images of the corresponding master frames (mf. 1 in green, mf. 2 in magenta, overlap in white). Small translations or rotations quickly lower the cross correlation, as in schematic in Fig. 1e, bottom.

Averaging all aligned master frames in a dataset produced a single “queen frame”. Using local maxima, we partitioned the queen frame into non-overlapping sectors, each representing the space that a single chromatophore can occupy (Fig. 1f). The number of chromatophore-pixels in each sector was then used to quantify the expansion state of the corresponding chromatophore in each image of the dataset (Fig. 1g). We thus obtained tens of thousands of parallel and simultaneous chromatophore-activity times series, characterizing skin patterns and their evolution (Methods).

Cuttlefish chromatophores come in different colours (Fig. 2a) usually classified in three to five groups9,18,19. To characterize chromatophore colour objectively, we measured their transmission spectra in freshly dissected skin (Fig. 2b, Methods). The distribution of spectra at 615 nm was bimodal, with a “dark” and a “light” cluster, plus intermediate colours ranging from orange to red. This could be explained in part by expansion state: local application of L-glutamate13, which induces chromatophore expansion, caused spectral changes towards lighter colours (Extended Data Fig. 3), consistent with previous descriptions20,21, and possibly explained by decreased local density and nano-structural features of the pigment granules22. Using techniques from analytical chemistry, we identified xanthommatin as a pigment in Sepia skin, and localized it exclusively to light chromatophores (Extended Data Fig. 3). Therefore, we can segregate Sepia officinalis chromatophores into two groups (light and dark) defined respectively by the presence or absence of xanthommatin. This confirms our initial classification based on transmission spectra (Fig. 2b).

Figure 2

Classifying chromatophores by colour

a, Isolated skin sample illustrating range of chromatophore colours. b, Transmission spectra of chromatophores including those shown in a. Lines correspond to ROIs of individual chromatophores, coloured as the average RGB value for that ROI. Right: Frequency distribution histogram for transmission at 615 nm. c, Small patch of averaged, aligned colour image (N=28,998 frames) showing dark (♢)/light () colour assignment. d, Normalized RGB colour distribution for pixels centered on N=9,199 chromatophores on one animal’s mantle Symbols as in c. e, Radial-averaged density of L and D chromatophores centered on D (left), and L (right) reference chromatophores. N=11,535 D; 14,802 L chromatophores, from 3 animals.

Extended Data Figure 3

Identification and localization of xanthommatin in light chromatophores

a, Chromatophores before (left) and after (right) local application of 40 µM glutamate. b, Transmission spectra of a population of chromatophores before (circles) and after (squares) glutamate application, projected onto 2 dimensions defined by human L and S cone action spectra (N=63 chromatophores). c, Direct infusion electrospray ionization Fourier-transform ion cyclotron resonance (ESI-FT-ICR) -mass spectrum of skin tissue extract showing high spectral complexity. d, High-performance liquid chromatography coupled to ultraviolet absorption spectroscopy and mass spectrometry (HPLC-UV-MS) chromatograms of skin tissue extract showing two main peaks with correlating UV (250nm) absorption (blue) and MS intensity (grey) comprised of eluting compounds with m/z 380.09 and m/z 424.08 (EIC traces, green). Experiments were replicated 5 times with similar results. e, Direct infusion ESI-FT-ICR mass spectra of skin tissue extract showing an overlay of the isolated precursor spectrum for decarboxylated xanthommatin (green, m/z 380.0876, theoretical: m/z 380.0877) and the fragment spectrum (red). Main fragments were assigned to putative structural losses of A (-NH3), B (-H2O, - NH3), C (-NH3, -HCOOH), D (-C2H3NO2) by accurate mass. f, Direct infusion ESI-FT-ICR mass spectra of skin tissue extract showing an overlay of the isolated precursor spectrum for xanthommatin (green, m/z 424.0774, theoretical: m/z 424.0775) and the fragment spectrum (red). Main fragments were assigned to putative structural losses of A (-NH3), B (-H2O, - NH3), C (-NH3, -HCOOH), D (-C2H3NO2) by accurate mass. g, Intensity distributions in laser desorption ionization Fourier-transform ion cyclotron resonance mass spectrometry (LDI-FT-ICR-MS) imaging and structures for putative xanthommatin derivatives (merged [M+H]+/[M+Na]+): decarboxylated, oxidized (m/z 380.0886; 402.0696), decarboxylated, reduced (m/z 382.1037; 404.0853), oxidized (m/z 424.0785; 446.0629) and reduced (m/z 426.0938; 448.0761). h, Intensity distributions of main xanthommatin and derivative fragments, corresponding to molecular species detected in ESI-FT-ICR fragmentation measurements. Experiments were performed on 12 tissue slices, producing similar results. i, Image of cryotome section of fresh-frozen sepia skin showing chromatophores. j, Spatial segmentation map of section in i, showing distinct clusters for light and dark chromatophores (orange vs black colours), and surrounding tissue (grey). k, Intensity distributions for xanthommatin derivatives (merged [M+H]+ and [M+Na]+) obtained from LDI-FT-ICR-MS imaging experiments (N=1 sample).

Consistent with results in vitro, chromatophore colour in vivo defined two modes, with partially overlapping dark and light (yellow to brown) clusters (Fig. 2c,d, Methods). The spatial arrangement of chromatophore colour was not random (Fig. 2e): we calculated the average local density of each colour class centred on chromatophores of a single colour (Fig. 2e; 32,740 chromatophores, 3 animals, Methods). On average, chromatophores of either class occupied a ~20 µm radius area. Beyond this, the density of opposite-colour chromatophores increased and dominated, peaking at ~55 µm. At ~100 µm, colour densities were inverted, indicating an alternation (on average) between light and dark chromatophores (see also Extended Data Fig. 4).

Extended Data Figure 4

Chromatophore-centered average densities

Two-dimensional density distributions for L and D chromatophores over an animal’s mantle (N = 9,199 chromatophores). The composite images show light chromatophore density in green and dark chromatophore density in magenta. For visualization, densities were linearly scaled together within an image. This preserves relative densities within each image but leads to slightly different colours across images.

Decomposing chromatophore control

To infer the potential structure of the control circuitry of the chromatophores, we examined their temporal co-variation during spontaneous fluctuations in vivo. This analysis is complicated by the facts that each chromatophore may be multiply innervated and individual MUs may overlap14, 23–28 (Extended Data Fig. 5).

Extended Data Figure 5

Identification of motor-units

a, Schematic showing three hypothetical, partially overlapping MUs (defined by motoneurons B, R and G), tracked over 3 epochs (i-iii) each characterized by different co-activation patterns (epoch i: R alone, ii: B + R and iii: B + G). Even though chromatophores 1 to 4 all belong to the same MU (R), their average pairwise correlation during these 3 epochs would differ due to the activity of the partially overlapping MUs B and G; identifying MUs using this metric would thus fail. This toy example indicates that the units of coordination during behaviour could be smaller than single anatomical MUs. (They could also be larger, if for example some motoneurons are always centrally coupled.) b, Single trials of in situ nerve minimal electrical stimulation experiments. Composite images (one per trial), green: 10 ms pre-stimulus, magenta: 200 ms post-stimulus. White: overlap. Threshold stimulation either leads to the expansion of a set of 3 chromatophores (marked with red circles, e.g., trial 1), or fails to activate any chromatophore (e.g., trial 6, 114 MUs determined with this method). c, Colour assignment of chromatophores in situ. Colour label was assigned based on a threshold on the red channel of RGB space (0.3). Chromatophores (dots) belonging to the same MU (as determined in a) are connected by lines, revealing the monochromaticity of MUs. N = 114 chromatophores. d, Dark MEs tend to be larger than light MEs. QQ plot showing quantiles of the dark vs. light ME size distribution. Line = identity. e, Tail of distribution of ME spread is heavier with D than L chromatophores. QQ plot showing quantiles of the dark vs light ME spread (calculated as in Fig. 3d). Line = identity.

To identify MUs directly, we first carried out minimal electrical stimulation of distal nerve branchlets innervating freshly dissected dorsal mantle skin and measured resulting chromatophore expansion (Methods). Putative MUs were small (2-10 chromatophores), usually clustered (median radius of 2.5 chromatophores or 247 µm) and monochromatic (Fig. 3a, Extended Data Fig. 5), consistent with observations in squid and octopus12,21,27,28. Light MUs were harder to stimulate electrically in isolation than dark ones, suggesting smaller axons.

Figure 3

Inferring chromatophore neural control from co-variation

a, Summary statistics of in vitro experiments. Top: fraction of dark chromatophores in a MU. Middle: Number of chromatophores in a MU. Bottom: Average distance to MU centroid for all chromatophores in that MU, normalized to the average nearest-neighbor distance over all chromatophores (N=295 chromatophores, 114 MUs). b, Identification of motor elements (MEs) in vivo. Size-over-time traces for 18 chromatophores over 2 chunks. Red: 9 chromatophores clustered as one ME. Blue: Nearest-neighbours (in physical space) to the chromatophores in red. Calibrations: x=5 s; y=1000 µm2. c, Average colour image showing position and colour of chromatophores in b (circles, colours). d, Summary statistics of in vivo experiments. Compare to a. e, Aligned colour image showing “average” pattern over a ~1 h-long dataset (237,826 frames). f, Left: Correlation-based hierarchical clustering of average ME time-courses for the dataset in f (N=695 MEs). Right: fraction of monochromatic clusters as function of correlation distance (N=1,896 MEs, 3 animals). g-i, Clusters at threshold levels in f (arrows, top to bottom) within frame in e. Same symbols throughout; colours denote cluster identities.

Returning to behaving animals, close observation revealed pronounced coordinated fluctuations of small localized groups of chromatophores, suggesting common drive (Supplementary Video 3). Seeking to extract these groups statistically, we factorized the chromatophore activity matrix using independent component analysis (ICA)29,30 and clustered through thresholding chromatophores contributing strongly to single independent components, allowing for multiple-cluster membership (Methods). We called these inferred clusters of chromatophores “motor elements” (MEs) to distinguish them from anatomically defined MUs (see above).

Using 57±10-min long datasets, we could extract hundreds of MEs over an animal’s mantle. While chromatophores within a ME tended to be highly correlated with each other, we often observed subsets within a ME that occasionally fluctuated independently of the others. We also regularly saw transient co-fluctuations with unclustered, otherwise weakly correlated, chromatophores (Fig. 3b).

MEs were mostly monochromatic (89% contained only light or dark) and their size distribution was heavy-tailed, with a median of 3 chromatophores (Fig. 3d, Extended Data Fig. 5). Note that their size distribution resembles that for presumed MUs (Fig. 3a, d)—with the tail likely representing groups of highly coordinated MUs, identified by ICA. Chromatophores within a ME were typically clustered physically, but were usually not nearest-neighbours, consistent with colour alternation (Fig. 2e) and monochromaticity. The distribution of spatial clustering was also heavy-tailed, with few MEs containing chromatophores spread over large areas (Fig. 3d).

We next tested whether the statistical relationships between MEs might reveal elements of higher-level control9. We averaged the time series of all chromatophores in single MEs to approximate their underlying neural drive (Methods). We then performed hierarchical clustering on the correlation of these average time series, clustering chromatophores composing MEs at different levels of granularity.

We illustrate the approach with a dataset in which the animal displayed distinct macroscopic pattern components (Fig. 3e). The first bifurcation (Fig. 3f, top arrow) divided chromatophores within the white square + posterior spots9 from all the others (Fig. 3g). Within these two super-clusters, many degrees of correlation (i.e., putative levels of common control) and monochromaticity could be identified. At lower levels in the hierarchy (Fig. 3f, middle arrow), the decomposition revealed medium-sized but still spatially structured elements, such as those that define the borders of the white square (Fig. 3h). Notably, some of these sub-pattern elements were at times observed to vary independently of their super-cluster, consistent with previous brain stimulation experiments indicating multiple levels of motor control31 (Supplementary Video 4). In turn, these elements could be decomposed into MEs (Fig. 3f, bottom arrow), often forming common or co-linear patterns indicative of precise innervation honouring the borders of macroscopic pattern elements (Fig. 3i).

Tracking pattern dynamics

Changes in an animal’s visual scene usually triggered rapid skin pattern changes. In the example in Fig. 4, a hand was moved above the cuttlefish (Fig. 4a), causing it to transition from dark to light (Fig. 4a, Supplementary Video 5). We examined this transition over several repeats by tracking the states of 17,305 chromatophores at 60 images/s (Fig. 4b-d). Projected into principal-component space for visualization, the data took the form of looping trajectories joining dark (i) and light (ii) states through a sequence of intermediate states (e.g., iii, Fig. 4b). Upon each stimulus, the animal not only generated the same target patterns (i and ii) but moved with chromatophore-level repeatability (Extended Data Fig. 6) through very similar, low-dimensional sequences of intermediate states (Fig. 4c-e: 85% of the variance is explained by 9 dimensions). This reliability of sequential chromatophore activation is remarkable because no physical constraints—such as those imposed by a moving limb for example—exist to prevent arbitrary and possibly more direct transitions. This suggests that neural constraints, probably linked to the putative control hierarchy inferred above (Fig. 3, Extended Data Fig. 7) and to internal connectivity, limit the paths along which pattern transitions can occur.

Figure 4

Tracking pattern changes at cellular resolution

a, Snapshots (i-iii) of an animal reacting to motion. Left: raw images. Right: corresponding segmented images: chromatophores shown as disks proportional to actual size (N=17,305 chromatophores, see Supplementary Video 5). b, Full sequence of skin patterns (17,305-D vectors of chromatophore sizes, unfiltered) over repeats of the behaviour, projected into space of the first 3 principal components (PC 1-3). Time in colour. c, Full area-over-time matrix of the behavioural sequence in b. Chromatophore areas normalized for visualization (0-1) and ordered by time-to-cross mean activity during first sequence (15-88 s). d, Correlation matrix of full, 17,305-D vectors of chromatophore sizes (Pearson’s). e, Cumulative variance explained by increasing numbers of PCs (dimensions). Black: sequence in b-d, Magenta: 37-min dataset, including more patterns. Arrows indicate x where y = 85%.

Extended Data Figure 6

Pattern-border precision at single-chromatophore level

Left: Three similar points along the pattern trajectories shown in Fig. 4b after chromatophore alignment. Right: expanded view of a pattern border. Note the remarkably similar expansion states of the chromatophores at each one of the three visits, and the rugged pattern borders at chromatophore scale, with interdigitation of expanded and contracted chromatophores, generating apparent noise. This apparent noise may be critical for natural realism.

Extended Data Figure 7

Linking statistical hierarchy of pattern elements to dynamics

a, Three exemplar intermediate-level clusters of MEs (threshold of 0.4 as in Fig. 3i, separate animal), overlaid on the average aligned colour image for the dataset (216,160 images). The clusters are mostly composed of chromatophores of a single colour: cluster 1 (red) is light; cluster 2 and 3 (green and blue) are dark. b, The dynamics of a 60-min. dataset, projected onto the first 3 principal components (48% variance explained, N=1437 chromatophores, 52,040 samples). A cluster activity direction can be defined in PC space by projecting the cluster identity vector (vector of length = number of chromatophores, with 1’s assigned to chromatophores in a cluster, 0’s otherwise), onto the PCs. The coloured lines show the cluster activity directions for the 3 clusters in a (line colour). Projecting the dataset onto these directions shows the expansion strength of the cluster at different times. The images corresponding to the times of lowest and highest strengths are shown to the left and right, respectively. c, Full distribution of expansion strengths, projecting all time points onto cluster activity directions. In this dataset, cluster 2 is often expanded, whereas clusters 1 and 3 are rarely expanded.

Tracking array development

Sepia officinalis continuously add new chromatophores as they grow, increasing from a few thousands in hatchlings to a few millions before death9. To track chromatophore insertion and development, we aligned multiple datasets recorded days apart (Methods, Supplementary Video 6). We observed that all chromatophores change colour in a systematic progression: all newly born chromatophores were pale yellow (Fig. 5a; Extended Data Fig. 8), consistent with observations in hatchlings34. In a 7-day-old animal, yellow chromatophores transitioned over ~two weeks to orange and later briefly to red. 18.7 ± 1.1 days after detection, each chromatophore turned dark and remained so throughout our observation period, possibly due to xanthommatin polymerization (Fig. 5b; Extended Data Fig. 8; 13 chromatophores, 25 days). Hence, the intermediate colours of chromatophores result from at least two causes: their expansion state (Extended Data Fig. 3) and their age (Fig. 5).

Figure 5

Tracking development of the chromatophore array

a, Aligned images of the same small skin patch over days, illustrating birth and colour changes. Crosshair zoomed-in below. “Day 1” applied post hoc to day preceding detection. Scale bars: ~50 µm (nonlinear alignment). b, Top: colour evolution over days for chromatophore in a, plotted in Hue-Value space. Bottom: same for N=13 chromatophores. After white-balancing, some black chromatophores appear bluish due to noise and are not shown (lower half of colour wheel). Blue line: median colour over all chromatophores at same estimated developmental time (N = 13). c, Theoretical relationship between new chromatophore insertion rate and colour maturation (transition to black) for L/D ratio measured in juvenile over development (blue), and for distribution of L/D ratios measured over animals 1-8 months old (green shows ± 1 s.d. from the mean, N = 7 animals, see Supplementary Information). Note that experimental measurements (red) fall precisely on theoretical curve (horizontal error bars: s.e.m., from cross-validation). d, Distributions of size (left) and standard deviation of size over time (right) for yellow, red (transitional) and dark chromatophores. Transitional chromatophores are significantly smaller and less variable over time than either yellow or dark ones (p=3.5e-92, 9.7e-107, Kruskal-Wallis followed by Tukey’s HSD, N=1,413 yellow, N=214 red, N=1,468 dark). e, New chromatophores arise in gaps in existing array. (i, ii): same skin patch aligned 11 days apart. Red dots centered on chromatophores detected on day 12 but absent on day 1. (iii): Grey-scale shows distance to nearest older chromatophore in ii. f, Summary of insertion location statistics. Dark: distribution of distances to nearest old chromatophore. Yellow: same distribution, conditioned on location of new chromatophore insertion. Blue line: Probability ratio of yellow to black distributions, showing ~monotonic increase at increasing distances to nearest old chromatophore. (N=11,527 old; 1,412 new chromatophores, 2 animals). NN: nearest neighbour.

Extended Data Figure 8

Chromatophores change colour from L to D as they age

A gallery of aligned patches of skin centered on the position of chromatophore insertion. Top: Juvenile animal, 7 days old on the first day of observation (D1). Left-most column shows skin pre-chromatophore-birth. Over ~19 days, chromatophores that first appear pale yellow darken progressively, transitioning to orange and red, before finally turning black. FOVs: from ~150x150 µm at D1 to 300x300 µm at D25. Bottom: Adult animal, 105 days old on D1. Rows show chromatophores undergoing a similar light-dark colour transition as in the juvenile above, but at a much slower rate. FOVs: ~200x200 µm (nonlinear alignment). Exemplars chosen from aligned skin patches containing ~100 chromatophores.

Chromatophores (i) do not disappear32 and (ii) their colour ratio (light/dark) is roughly constant (1.06±0.19:1 in seven 8-252 day-old animals) 33. However, (iii) the time over which chromatophores turn from light to dark increased from ~19 days (above) to ~97 days in a 105-day-old animal (96.6±9.3 days, s.e.m., Methods). Likewise, whereas light chromatophores are produced daily as a fixed fraction of all existing chromatophores, (iv) the rate of chromatophore addition dropped from 4.1% in a 7-day-old animal to 0.6% in a 105-day-old animal. In Supplementary Information, we provide a formal derivation of an expression linking these observed quantities (ii-iv). Figure 5c shows the theoretical interdependence of two of these quantities (iii vs. iv) given (ii): the lifetime of the light state and insertion rate measured experimentally fall precisely on this curve, suggesting that these two properties are balanced to maintain a near constant colour ratio across life span.

The monochromaticity of MUs (Fig. 3a) could result from the fact that new motoneurons9 innervate only newly born (light) chromatophores21. This hypothesis, however, introduces a conundrum in that each animal should keep track of the age of each MU to know its colour, an unlikely feat. This problem could be solved, however, if the chromatophores were re-innervated as they change colour, replacing “light” with “dark” motoneurons. Consistent with this, we observed that the average diameter and size-variance over time of the red chromatophores were smaller than with the light and dark ones (Fig. 5d, Extended Data Fig. 9), suggesting that MUs undergo re-innervation as they near transition to dark.

Extended Data Figure 9

Development of the chromatophore array

a, Flowchart depicting the spatial-growth-model algorithm and highlighting the involvement of model parameters (Methods). b, Boxplots of nearest-neighbor (NN) distances between young (<6 days old) and older chromatophores. Young chromatophores are significantly closer to both older light (>12d) and dark chromatophores than to other young or middle aged (6-12d) light chromatophores. (p<0.0001, Kruskal-Wallis followed by Tukey’s HSD, Nchromatophores = 522 <6d, 541 6-12d, 1550 >12d, 1910 dark, 1 animal). Distances calculated on single image, ages estimated by finding the day of chromatophore birth on aligned developmental datasets (Methods). c, Distributions of size for yellow, red (transitional, T) and dark (D) chromatophores, annotated manually (validation of analysis in Fig. 5d). Transitional chromatophores are significantly smaller than either yellow (Y) or dark ones (TvY, p= 1.0e-7, TvD, p= 6.3e-4, N = 70 Y, 16 T, 84 D, two-tailed Wilcoxon Rank Sum tests, 1 animal). Boxplots- central line: median, box limits: quartiles, whiskers: ±2.7 s.d. d, Generation of the inhibitory surround used in the skin growth model (Fig 6b). Blue: Empirical radially averaged chromatophore centered density, inverted and normalized 0:1. Red: Logistic function fit to the blue density, as in Fig. 6b. e-h, Manipulating single parameters of the skin growth model suggests the mechanisms underlying colour interdigitation. e, Difference between peak D-triggered D density and D-triggered (light chromatophore) L density, as a function of model skin growth rate. Points are from the average of 3 model runs. Line is linear fit. ANOVA F statistic vs. constant model=96.6, p=0.000186. f, Difference between peaks of radially averaged D-triggered D density and D-triggered L density, as function of age at which chromatophores transition from L to D. Points are from a single model run, where the colour class was changed according to chromatophore age. Line is linear fit. F-test for linear regression: F=152, p=5.26e-6. g, Difference between first peak (first zero-crossing of derivative of radially averaged density) in the radially averaged L-triggered L density and L-triggered D density, as a function of r, the rate at which the inhibitory surround changes with chromatophore age. Points are from the average of 3 runs of the model. Line is linear fit. F-test for linear regression: F=21.9, p=0.00226. h, Colour interdigitation is robust to stop-criterion used to define end of “day” (parameter 5, Methods). Magenta: D-triggered D density - D-triggered L density (as in e, f). Black: L-triggered L density minus L-triggered D density peaks. Lines are linear fits. F-test for linear regression: F=0.0206, 6.57, p=0.889, 0.0334. Points in e-h are from the average of 3 model runs.

We next examined the geometry of new-chromatophore insertions. The example in Fig. 5e i, ii shows the same aligned patch of skin at 11 days interval. The positions of the young (yellow) chromatophores at day 12 (red dots) retrospectively identified the zones of their future insertion as positions far from already born chromatophores (Fig. 5e iii, 5f). This arrangement suggested a simple model for the development of the chromatophore array, based on the regulated production of a hypothetical inhibitory signal by each chromatophore19,35, which we evaluated using computer simulations.

Simple rules can explain spatial layout

Our simulation ran on discrete steps (days) and was initiated by the random insertion of light chromatophores (L) in a patch of bare skin, constrained by a chromatophore-centered inhibitory surround (Fig. 6a, b). Once filled with chromatophores, the skin patch “grew” isotropically by a fixed proportion, followed by the next “day” of chromatophore insertion. When chromatophores reached 19 days (above) they switched from L to dark (D) (Fig. 6a). The inhibitory surround was described by a sigmoidal function derived from empirical measurements. To match the experimentally measured differences in spacing between newer and older chromatophores (Fig. 2e, Extended Data Fig. 9), we allowed the size of the inhibitory surround to change as chromatophores age. We fixed the shape of the surround, and fit its initial size, s, as well as the rate of size change with age, r to empirical measurements (Methods).

Figure 6

Simple rules explain the spatial layout of chromatophores

a, Schematic of model. Every day, new light chromatophores are inserted into skin, chromatophores move apart as the skin grows, and chromatophores turn from light to dark when they reach maturation age (19 days). b, Left: disc of surround-inhibition centered on each chromatophore. Right: radius of inhibitory surround decreases with age at rate r from an initial birth size s. c, Percentage of light chromatophores as function of simulated day (50 runs). Red line: experimental ratio. d, Real vs. simulated chromatophore position and colour assignment for small skin patches. e, Summary of insertion location statistics as in Fig. 5f, but for simulated developing skin. “New” insertions are those on the last “day” of simulation. (Four of 15,494 new insertions at 69-73 µm are not shown for clarity.) f, Distribution of Pearson’s correlation coefficients between a random patch and all locations of simulated skin. Inset: cross-correlation function centered on the correct matching location, corr. = 1. g, Radially averaged density of light and dark chromatophores from the centre of dark (left), and light (right) chromatophores. Stippled lines: experimentally measured densities (N = 4,095 dark; 5,104 light chromatophores). Thin lines: 50 independent runs of the model.

Consistent with our analytical result (Fig. 5c), this simple model converged to the observed percentage of light chromatophores (model=0.55±0.01, 50 simulations; data=0.55, 5104 L/ 4095 D, 1 animal; Fig. 6c), provided the skin growth rate was set to allow a realistic rate of chromatophore insertion (model=4.23±0.01%; data= 4.1%/day). It produced realistic spatial patterns of L and D (Fig. 6d), new chromatophore insertion locations (Fig. 6e vs. Fig. 5e) and chromatophore density (mean density chr./µm2: model=2.52e-4±0.1e-4; data=2.44e-4). Local patches of simulated skin had unique spatial layouts (“fingerprints”) of the type that we exploited for image registration (Fig. 6f vs. Fig. 1e). Our model was able to produce realistic local chromatophore-centred densities, featuring the experimentally observed inter-digitation of colour-specific modes (Fig. 6g vs. Fig. 2e, Extended Data Fig. 9). Notably, by varying r, we could generate other known chromatophore distribution patterns such as the discoid units observed in some squid species36 (Extended Data Fig. 10). This simple model may thus apply more generally to cephalopod skin patterning.

Extended Data Figure 10

Exploration of developmental-model parameters reveals species-specific patterns

Changing model parameters (see text and Methods) can lead to the characteristic rings observed in some squid species, with single light chromatophores at the center and a radial centrifugal darkening gradient. Top: skin of common squid, Loligo vulgaris (image: Robert Siegel). Bottom, simulation of development using a profile of change of the inhibitory disc centered on each chromatophore r different from that used in Fig. 6b for Sepia officinalis.

Discussion

We developed a strategy to track tens of thousands of individual chromatophores in freely behaving cephalopods, enabling studies of behaviour and development at cellular resolution. Our results open a path towards addressing many important biological questions. A first concerns visual perception. Cephalopod camouflage is unique in revealing a high-dimensional neural readout of an animal’s visual texture perception. Identifying the primitives of cephalopod camouflage might indeed reveal fundamental features of texture generation but also of vertebrate texture perception, since the former (in cephalopods) likely evolved in response to the latter (in their vertebrate predators). A second concerns the development of methods to analyse very large neural datasets in the context of naturalistic behaviour37. Because chromatophore data can be assigned unambiguously to identified elements that lie at the same level of a neural hierarchy (here exclusively motoneurons), their analysis does not suffer from assumptions about their identities and positions in structured or recurrent circuits, as may happen with brain neural imaging. A third concerns morphogenesis and development. Our data suggest that simple local rules can explain the structure of a continuously growing chromatophore array. They thus lead directly to clear questions about mechanisms and about their similarity with ones known from other systems38–40. Fourth, our results indicate that very complex behaviours can be described quantitatively at cellular resolution and in species that may reveal much about shared constraints on brain evolution41. This system is therefore particularly well suited for studying the relationship between neural and behavioural dynamics, a central and general problem in neuroscience.

Methods

Experimental animals

Animal experimentation in this study was performed according to German animal welfare law (paragraph 11, sentence 1, #1, German animal welfare law to house and breed cephalopods for scientific purposes). European cuttlefish Sepia officinalis were hatched from eggs collected in the English Channel and reared in a seawater system, at 20°C. The closed system contains 4,000 litres of artificial seawater (ASW, Instant Ocean Inc.) with a salinity of 33 ‰ and pH of 8-8.5. Water quality was tested weekly and adjusted as necessary. Trace elements and amino acids were supplied weekly. Marine LED-lights above each tank provided a 12/12-hr light/dark cycle with gradual on- and off-sets at 7am and 7pm. The animals were fed live food (either Hemimysis spp. or small Palaemonetes spp.) ad libitum 3 times per day. Experimental animals were selected for healthy appearance and calm behaviour. The animals were housed together in 120-liter glass tanks with a constant water through-flow resulting in five complete water exchanges per hour. Enrichment consisted of natural fine-grained sand substrate and seaweed (Caulerpa prolifera).

In vivo behavioural data acquisition

For in vivo behavioural experiments, 6 animals (1 to 60 days post-hatching, ~6-50 mm in mantle length) were placed in a capped filming chamber (150 mm × 95 mm × h75 mm, or 240 mm × 170 mm × h50 mm) filled with seawater. A single filming session typically lasted between 10 and 90 minutes per day and per animal. Our filming procedures induced no pain, suffering, distress or harm to the animal. Naturalistic textures with normalized power spectra (Normalized Brodatz Texture database) and artificial patterns generated in MATLAB and Paint were displayed on the floor of the tank using an E-Ink display. Filming was performed at 59.94 frames per second (fps) in 4K full-frame (4096 × 2160) using the Sony PMW-F55 camera in the Sony RAW format. Resolution was 40.8 ± 32.4 µm2/pixel. The camera was mounted on a motorized X-Y translation stage and its position adjusted with a joystick to keep the animal in view. Acquired data were colour matched in DaVinci Resolve Studio 12.5 (Black Magic Design). Movies were compressed offline to the H.264 format using the x264 encoder provided by FFmpeg-2.8.6, with the compression preset “faster” and the constant rate factor of 16, without chroma subsampling. These compressed movies were used for all subsequent in vivo data analysis. Note that the analyses and results presented in this publication do not depend on the exact statistics of the images or patterns shown to the animals. They depend singly on our ability to detect changes and correlations between patterns produced by single animals over time, at sub-chromatophore resolution.

In vitro electrophysiology

For in vitro experiments, animals (19 to 40 weeks old, 115 ± 28 mm in mantle length, 10 in total) were euthanized according to well-established best practice protocols42: animals were deeply anaesthetized first in isotonic 3% ethanol in ASW and then in 5% ethanol in ASW or using 3.5% (w/v) MgCl2 in ASW. Superficial skin samples were then removed gently from the dorsal mantle, peeling away the superficial skin layers from the underlying body musculature, gently cutting nerves and connective tissue with iridectomy scissors and taking care that the chromatophores were not overly stretched or damaged. These skin samples were placed in cold ASW inside a transparent observation chamber, pinned at their edges, superficial surface down, stretched gently so as to eliminate wrinkles and left to recover. The chamber was placed on the translation stage of an inverted microscope and the chromatophore array observed with 1.25× or 2.5× objectives (126.75 or 291.78 px/mm). A fine suction electrode operated with a micromanipulator was placed on the cut end of a nerve and for electrical stimulation (pulse duration: 100 µs) at pulse rates of 0.5 Hz and at threshold intensity using a pulse generator (± 5.57 to ± 2774 µA dialled gradually, A.M.P.I. Master-8-cp) and a constant-current stimulus isolator (World Precision Instruments A360). Colour images of the chromatophores were acquired with a CCD/CMOS camera (Basler acA1920-155uc) at 30 fps. The stimulus trace and camera exposure times were recorded using a digitizer (Axon Digidata 1440A). Synchronization and analysis were conducted off line.

Transmission spectra

Transmission spectra were recorded from fresh skin samples (extracted as above). Samples were mounted on standard microscope slides with #1.5 coverslips in ASW. Transmission spectra were recorded in 32 channels on a Zeiss LSM 880 Examiner confocal microscope (10x, NA 0.45, water immersion objective) using the “lambda mode” spectral detection. This mode is usually intended for fluorescence detection, but the halogen lamp for transmitted-light mode can be turned on with a service macro. Images were thus recorded with a scanned point detector but with widefield illumination. We manually drew regions of interest around chromatophores. Raw spectra were normalized with respect to a nearby region in the same field of view that did not contain any chromatophore. To measure the effect of expansion state on transmission spectra, ASW was replaced by a glutamate solution (40µM) in ASW. Images were acquired before and 6-8min after glutamate application.

Mass Spectrometry

Excised skin tissue samples were homogenized in 1:1 methanol:water (v/v) supplemented with 1% trifluoroacetic acid, sonicated for 10 min in an ultrasonic bath and placed for 2 h on a rotary shaker. The samples were then centrifuged for 10 min, the supernatant removed, filtered through a 0.2 µm syringe filter and evaporated to dryness. For mass spectrometry, the dried extracts were re-suspended in 95:5 water:acetonitrile (v/v) supplemented with 0.1% formic acid.

High-performance liquid chromatography coupled to ultraviolet (UV) -absorption and mass spectrometric detection (HPLC-UV-MS) experiments were carried out on an Ultimate 3000 RSLC system (Dionex, Germering, Germany) equipped with a CSH C18 column (Charged Surface Hybrid, 2.1 x 100 mm, 1.7 µm particle size, Waters, Winslow, UK) and variable wavelength detector (VWD) set to 250 nm, coupled to an Impact II mass spectrometer (Bruker Daltonik, Bremen, Germany). Separation was carried out using water (A) and acetonitrile (B), both supplemented with 0.1% FA, as mobile phases with a flow rate of 300 µl/min at 40°C. After 2 min of equilibration with 2% B, a linear gradient was ramped from 2% to 95% B in 30 min followed by 5 min wash (95% B) and 3 min equilibration (2% B) steps. The mass spectrometer was operated in positive-ion mode with a mass range of m/z 50 to m/z 1000. Processing and data analysis were performed manually using DataAnalysis 4.4 (build 200.55.2969, Bruker Daltonik).

For direct infusion experiments, extracts were diluted 1:100 and infused at 120 µl/h into a 7T SolariX XR mass spectrometer (Bruker Daltonik). Spectra were recorded in positive-ion mode of m/z 107.5 to m/z 2000. For exact mass determination and fragmentation experiments, precursor ions were isolated using the quadrupole, inspected for contaminating ions and then subjected to CID in the collision cell. Spectra were analysed with DataAnalysis 4.4 (build 200.55.2969, Bruker Daltonik).

For MSI experiments, excised skin tissue samples were stretched and pinned onto frozen gelatine blocks, snap-frozen in isopentane and sectioned to 12 µm using a CM 3050s cryotome (Leica Biosystems, Nussloch, Germany). The slices were carefully transferred onto conductive ITO-coated glass slides (Bruker Daltonik), thaw-mounted and dried in a vacuum desiccator prior to taking optical slide scans with an OpticLab H850 histology slide scanner (Plustek, Ahrensburg, Germany). Samples were screened using a Rapiflex TOF/TOF mass spectrometer (Bruker Daltonik) operated in positive- and negative-ion modes, using a mass range of m/z 100 to m/z 2000. Ultrahigh-resolution mass spectra were acquired on a 7T SolariX XR mass spectrometer (Bruker Daltonik) in positive-ion mode in a mass range from m/z 107.5 to m/z 2000 using a 20x20 µm pixel grid. The laser was operated at 500 Hz with 100 shots/pixel and focus set to minimum. Imaging data were acquired and pre-processed using flexImaging 5.0 (build 5.0.78.0_1031_152, Bruker Daltonik) and further analyzed and visualized using SCiLS Lab 2016b (build 4.01.8758, SCiLS, Bremen, Germany). Individual images were adjusted to the same intensity scale and weak spatial denoising applied for merged compounds. Spatial segmentation was performed with weak spatial denoising and a bisecting k-means algorithm based on the correlation distance of individual spectra. The relationship between colour and xanthommatin concentration was examined by manually clustering partitions of the k-means algorithm corresponding to yellow and red/brown chromatophores. We then took the average xanthommatin concentration in each cluster, summed over all derivatives.

Summary of the image processing and tracking pipeline

The major steps of the processing pipeline were:

Chunking: Identify episodes of video (‘chunks’) with cuttlefish in focus.

Segmentation: Label pixels as chromatophore or background on individual frames.

Registration: Alignment across frames within chunks to correct slight nonlinear body distortions (over seconds).

Stitching: Alignment across chunks (seconds to hours)

Chromatophore identification and size tracking.

Colour assignment.

Stitching across days.

Chunking of in-focus frames

In vivo behavioural datasets consisted of series of frames where the animal was in view and in focus, separated by frames where the animal was out of view, out of focus or blurred due to fast motion. We first identified in-focus frames using a simple focus statistic (sum of a difference-of-Gaussians filter size 15x15 pixels, σ = 1.5 and 2 pixels) to each image. The standard deviation of the filter was selected to match the mean size of chromatophores. Our statistics therefore indicated whether chromatophores were present and clear in an image. Continuous sequences of images were then selected semi-automatically, based on the amplitude of the focus statistic and its variability over images and time. We called “chunks” such continuous in-focus image sequences obtained within a single filming session.

Chromatophore Segmentation

We segmented chromatophores from the background using a supervised learning approach. For training and validation, images (1 each) of 256 x 256 pixels containing a representative sampling of chromatophore sizes and colours were manually annotated pixel-wise as belonging to a chromatophore or background skin. Annotations were performed by six individuals and inconsistencies were removed using majority vote. We then fitted a random-forest model16 to this annotation. Our model classified pixels as chromatophore or background based on feature vectors calculated from the output of 8 difference-of-Gaussian filters (σ = 0.81.4x, x = 1:8) per RGB colour channel. Filter sizes were chosen to cover the range of observed chromatophore sizes. We determined the random forest parameters by hyper-parameter optimization43 (number of trees in [1, 32], maximal depth in [1, 32], minimal data-size for split in [1, 11], minimal data-size for leaves [1,21], splitting criterion either by Gini impurity or information gain and enabling or disabling bootstrap aggregation). 1,000 models were fitted using a 4-fold cross-validation and the best model was identified by the F1-score. To this end we used a model with 8 trees with a depth of 8 and an entropy-based splitting criterion on 5 randomly selected features. We assessed model performance by comparing classifier performance against a second, manually annotated image (Extended Data Fig. 1).

Registration: alignment of images within a chunk

Animal movement (e.g., breathing) and skin distortions caused the pixel location of chromatophores to change over frames. Our high frame rate combined with the definition of chunking meant that differences in chromatophore locations (both affine and non-affine deformations) between successive frames within a chunk could be assumed to be small. We could therefore use image registration methods for small-displacement optic flow. We used the Lukas-Kanade optical flow algorithm17 to track points centred on a random subset of ~300 round chromatophores. Round chromatophores were found by placing a threshold on circularity of chromatophores detected in the first frame of the chunk. These chromatophores were selected to minimize runtime. The full-frame optical flow was interpolated from these tracking points using a moving least-squares algorithm44. We chose a smoothness parameter (alpha = 3) for interpolation to remove skin distortions and large movement, but not the fine scale movement of individual chromatophores.

Stitching averaged aligned images over chunks

Once all the images within a chunk were aligned, we averaged over the binarized images, generating a “master frame”. The value of each pixel in a master frame thus represents the fraction of frames within the corresponding chunk in which that pixel was labelled as belonging to a chromatophore. Because chromatophore size can vary over frames during a single chunk, the typical profile of a chromatophore in a master frame is a radial gradient. After obtaining one master frame per chunk, we developed a method to register all the chunks of a filming session into a common reference frame. We call this process “stitching”. Individual chunks were, by definition, separated from each other by out-of-focus epochs, i.e., ones in which the cuttlefish often changed position, angles in x, y and z, body shape, and chromatophore pattern. Stitching thus required aligning and morphing chunks into the same reference frame. For every master frame in a dataset, we first defined a mask outlining the cuttlefish by applying a difference-of-Gaussian filter to the image and thresholding the result. We then mapped all master frames into each other’s reference frames.

To stitch together 2 master frames A and B, we first performed a coarse rigid body transform mapping A into the reference frame of B by fitting an ellipse around the cuttlefish mask in both frames. This created B’, i.e., B mapped into the reference frame of A through the inverse of this mapping. Next, we defined a grid of points 256x256 pixels apart over the cuttlefish mask of A. We attempted to find each point of this grid in B’ by correlating patches of 64x64 pixels centred on the grid points in A with regions in B’. We sampled a range of translations (±256 px in 2-px steps) and rotations (±20° in 2° steps) around the pixel location of each grid point to find the pixel of highest correlation value. In general not all grid points could be correctly mapped; outliers were removed using the RANSAC algorithm45 under an affine model. A new map was constructed from the remaining mapped points using moving least squares interpolation44.

By applying the inverse of this mapping to B’ we produced B”, a more refined mapping. Fine alignment was performed by repeating this process on a finer grid. A new grid of points 32x32 pixels apart was defined on the cuttlefish mask of A. We then attempted to find each of the points in B” with the highest local cross-correlation to 64x64 pixel patches centered on each grid point. We then interpolated between these points using moving least squares to produce a full map, B”’. Combining these three maps resulted in a single non-affine mapping from B into the reference frame of A.

This stitching algorithm was used to map every master frame in a dataset into the reference frame of every other master frame. We could quantify the accuracy of this non-symmetric mapping over the cuttlefish mantle by calculating the reprojection error: all points in the cuttlefish mask of master frame A were mapped into the reference frame of master frame B (using the A to B map) and then mapped back into the reference frame of A (using the B to A map). The reprojection error was defined as the Euclidian distance between the original and remapped points. A point was considered well mapped if it reprojected to within 3 pixels (20±6 μm) of its original location. By taking the fraction of well-mapped points in every master frame we produced a matrix quantifying how well every master frame mapped into every master frame over the cuttlefish mantle. The column of this matrix with the highest well-mapped fraction identified the master frame into which all others mapped best. We used this as common reference frame for the dataset (see below). Poorly mapped chunks of data, defined as master frames that had less than 50% well-mapped points within the cuttlefish mask were excluded from subsequent analysis. Stitching failures were usually the result of poorly registered chunks, resulting in blurry master frames. These failures, in turn, were often due to temporary loss of focus through cuttlefish moving in and out of the focus range of our optical system.

Chromatophore definition

As done for within-chunk alignment, we used the maps generated by our stitching algorithm to project all master frames into the dataset’s common reference frame. The resulting average frame (“queen frame”) represented approximately the probability of each pixel being labelled as belonging to a chromatophore throughout all in-focus frames over the entire filming session, excluding poorly mapped chunks of data (see above section). Chromatophores were detected by finding local maxima using a 3x3 footprint kernel. Applying the watershed transform to the queen frame with detected chromatophores as markers divided the image into basins, each defining a region surrounding a single chromatophore. The watershed transform also split groups of merged chromatophores, relying on gradients in the queen frame created by chromatophore size changes over the dataset. A conservative mask defining the region of the cuttlefish mantle that was in focus was drawn as a convex hull of manually selected points, and all subsequent analysis was performed on chromatophore basins within this masked region.

Chromatophore-size tracking

With the chromatophore basins so defined, we could track the expansion state of each chromatophore over time. Each segmented image (see chromatophore segmentation, above) was mapped into the reference frame of the first image in a chunk (see alignment of images within a chunk, above), and then mapped again into the common reference frame of the dataset (see stitching average aligned images over chunks, above). The number of pixels classified as belonging to chromatophores in the segmented image was counted in every watershed basin of the queen frame. This pixel number, multiplied by an experiment-specific µm2/pixel2 calibration, defined the area of each chromatophore in each frame. This calculation was repeated with all images in a dataset, producing parallel time-series measurements of size over time for all segmented chromatophores. While chromatophore size generally varies with mantle position as the animal adopts different skin patterns, a uniformly dark pattern, as seen in Fig 4a, reveals no strong correlation between chromatophore size and anterior-posterior/medial-lateral position (r = 0.06, -0.03).

Colour assignment

Determining the colour of chromatophores is difficult to accomplish accurately on single images in vivo due to camera pixel noise, variability in lighting, and the expansion-state dependency of chromatophore colour. We therefore analysed chromatophore colours by mapping all images of a dataset into the dataset’s common reference frame, producing an average colour image. We first constructed a feature space where chromatophores could be accurately colour-labelled independently of our segmentation algorithm. We high-pass filtered the image and then performed independent component analysis in color space29,30. After projecting the image onto the two largest independent components and thresholding each projection separately, we took the maximum value over projections. This image was smoothed with a Gaussian filter (s.d. 1 pixel). We then applied a watershed algorithm46 to identify chromatophore regions. Chromatophore centers were defined as the weighted centroid of each region. Visible chromatophores that were not detected automatically, typically smaller red and yellow chromatophores, were identified manually. The average RGB value of a 3x3 pixel region of the average colour image, centred on each chromatophore’s location, defined chromatophore colour. These colours were clustered into 2 classes by fitting a Gaussian mixture model. We observed no strong correlation between the anterior-posterior or medial-lateral position on the mantle and chromatophore colour label (r=-0.0047/-0.0029). We defined the colour of tracked chromatophores by performing a nearest-neighbour matching between chromatophore centres defined using our tracking pipeline and centres defined on the average colour image. For the ME inference experiments, 93% of dark and 82% of light chromatophores detected in the average colour image (N = 39,948) could be linked to a tracked chromatophore (N = 35,062) located within 3 pixels (15.1 ± 5.9 µm, from 3 animals).

Stitching over filming sessions and detecting new chromatophores

We filmed 2 additional animals over periods of several weeks, i.e., periods over which the animals underwent considerable growth. To track chromatophores over days and weeks, we used a modified version of our stitching algorithm (see stitching average aligned images over chunks, above). To model cuttlefish growth, we used a similarity transformation rather than a rigid body transformation in the initial coarse alignment of cuttlefish masks. In the subsequent two correlation-matching steps, we use a larger search space in scales, rotations and translations. We used the resulting map to warp one dataset’s queen frame into the reference frame of another’s. After this we linked chromatophores over days through a nearest-neighbour matching. A convex hull containing the intersection of both dataset’s chromatophore basins was first calculated, and matching was only performed within this region. New chromatophores were defined as chromatophores from the later dataset that were located within the convex hull intersection and were not matched with a chromatophore in the earlier dataset. 75±14% of the pixels in the convex hull that contained mapped chromatophores had reprojection errors below 50 µm (6 datasets mapped from 3 animals), allowing for unambiguous matching through manual inspection and modification of matches using a custom GUI.

Pipeline implementation

Our alignment and tracking pipeline was implemented on two computing clusters: The Draco supercomputer at the Max Planck Computing and Data Facility (Garching, Germany) where 16 jobs were processed in parallel on 1-2 nodes with 32 cores @2.3GHz (128 GB RAM/node); and the FIAS computing cluster where 3 jobs were processed in parallel on 2-4 nodes with 32-64 cores/node @2.6 GHz (64 GB RAM/node). Data management between local storage and compute nodes was managed by a Bash script determining the sending and receiving of data and configuration files, and starting the pipeline on compute nodes. On each compute node, the pipeline computations were managed by a second Bash script which inserted all pipeline steps into a SLURM47 queue. Parallel computation of steps was handled by the SLURM controller. Parallelization per step was achieved by spawning programs using MPI48, and distributing computations across program instances. For steps where random access across video frames could not be implemented, e.g., reading a video, MPI spawned programs followed a one-producer/multiple-consumers pattern. For registration, parallelization was achieved by distributing the optical flow and moving least-squares algorithms across program instances. Threading was done using third-party libraries. The pipeline was written for GNU/Linux operating systems in the Python and Cython programming language, relying on Scipy49, scikit-learn50, scikit-image51, PyQt, and OpenCV-Python. Further data analysis was performed in Python and Matlab. The results of all pipeline steps were stored using the HDF5 format. The pipeline was constructed so that each step read one file and output another. In total, we achieved an overall speed of ~1-1.5 frames/s (corresponding to ~2.4 - 3.6 MB/s of a compressed video).

Chromatophore-triggered density plots

We constructed images composed of the locations of certain chromatophores (light, dark, etc.) and then averaged regions of these images centred on the location of chromatophore classes of interest. The resulting chromatophore-triggered average image was linearly interpolated to 1 px/µm2 and then smoothed with a Gaussian filter (15 µm size, 2-3 µm s.d.) We then computed radial averages.

Automated extraction of motor units (in vitro experiments)

Consecutive frames were first aligned using optic flow (as above) to correct for spontaneous skin movements. For experiments containing >7 consecutive stimulus trials, the video frames nearest 10 ms pre- and 200 ms post-stimulus were inspected and annotated manually for expansion events using image subtraction. Chromatophore position and colour was determined from the pre-stimulus frame of the first expanding trial. Colour classification was made using a threshold on the red channel of white balanced, contrast stretched RGB space. MUs were identified by coincident responses and failures of a set of multiple chromatophores (>1) across trials, allowing for an individual failure rate of up to 25%. We estimated the average spontaneous expansion probability (0.0293) by examining the activity of 81 chromatophores from 4 animals at times without stimulation. We then could estimate the chance probability of an observed sequence of responses and failures as where m is the number of chromatophores in a putative MU, e is the number of expansion trials, and n is the number of total trials. A threshold of p=0.05 was placed on observed sequences for inclusion as a MU. MUs along the edge (mean + 1 s.d. of the average chromatophore nearest-neighbor distance) of the field of view were excluded from the analysis to prevent underestimation of MU size.

Inference of Motor Elements from in vivo imaging data

Our choice of statistical model for motor unit inference was motivated by the desire to capture the potentially overlapping innervation of motor units while excluding sets of chromatophores that are more transiently coordinated. The non-Gaussian nature of the area fluctuations, at times suggestive of sparse motor neuron activity suggested that independent component analysis29 (ICA) would be well suited.

Chromatotophore-area time series were symmetrically low-pass filtered to 4 Hz using a 3-pole Butterworth filter. They were then down-sampled 4-8 fold. We performed “spatial” ICA on the resulting matrix, using the Fast ICA algorithm30. This algorithm iteratively estimates S = W X, where X is the centred, whitened, area-traces chromatophore matrix, W is the unmixing matrix, and S is the component chromatophore matrix of independent components. We used the algorithm to estimate C independent components, where C is the number of dimensions explaining 99.5% of a dataset’s variance.

To estimate statistically small sets of chromatophores receiving common drive (‘motor elements’ or MEs, see main text), we subsequently clustered the small subset of chromatophores with high values on single independent components (ICs) of the matrix S. Since motor-unit membership is binary (a chromatophore either is innervated by a motor neuron or is not), we thresholded the ICs to extract these highly contributing chromatophores and examine their properties. We found that the highest contributing chromatophores most often clustered spatially in single modes, with chromatophores contributing less located further away. We chose our threshold for ME inclusion such that the median of the chromatophore spatial distribution matched approximately that measured in vitro. We assigned a sign to each IC as the sign of the maximum value (chromatophore) of that IC. Values higher than 8 s.d. above the mean value of positive ICs, or lower than 8 s.d. below the mean value of negative ICs were clustered to form a ME. We then visually inspected the MEs to check for colour classification errors and to remove groups that did not contain well-segmented chromatophores due to errors in watershedding.

Inference of putative motor control hierarchy

We first averaged the filtered, down-sampled area time series (as in Motor Element inference above) for all chromatophores assigned to a ME (ignoring the sign or weight of its associated IC). This procedure was motivated by the known underlying biology: it attempted to approximate the common motoneuron drive that caused the chromatophores to be clustered into a ME. Note that the precision of this approximation depends on several factors, including the multiple and partially overlapping innervation of chromatophores, the difficulties of inferring motor units (as described above), and the fact that the relationship between chromatophore size and neural drive is likely sigmoid, and thus linear only in a limited range. We performed agglomerative hierarchical clustering of these ME time series, using the correlation coefficient as a distance metric and complete linkage. To segment monochromatic clusters at different levels of the hierarchy, we measured the fraction of clusters composed of MEs that contained only light or only dark chromatophores.

Chromatophore colour changes over development

For precise characterization of chromatophore colours (Fig 5a,b), we took images of 2 cuttlefish over days using an 18 M-pixel camera (Canon, 550D) at 10-18x magnification. Recognizable landmarks (e.g. papillae, mantle edges) were used to return to the same area of skin repeatedly. We aligned skin patches using TrakEM2 (ImageJ plugin). Images of chromatophores were white-balanced using a nearby patch of skin that did not contain any chromatophore as reference. The colour of individual chromatophores was then determined in hand-drawn regions of interest (ROIs). To visualize the colour change in a condensed representation in colour space, the three colour channels (red, green, and blue) were averaged over all pixels within the ROI. Colours were then converted from RGB into HSV colour space, which assigns brightness and hue to different axes. To determine an average trajectory in colour space, the chromatophores were temporally aligned to the transition state via thresholding on the red colour channel.

Chromatophore sizes over development

To check for potential size and variability differences of the transition state, we aligned a dataset separated by 2 days using our pipeline and identified transitioning chromatophores as those that were classified as light on the earlier day and dark on the later one. Yellow chromatophores were defined as light chromatophores that were not transitioning. Size and variability were estimated using filtered data (as in ME inference). For validation independent of our tracking pipeline (Extended Data Fig. 9c), we aligned images using TrakEM2, and manually grouped individual chromatophores into three categories (yellow/orange, reddish-brown, black) without temporal context. We then incorporated developmental information by retaining only those yellow/orange chromatophores that were observed again at a later time point as yellow/orange, i.e. those data, that were not close to the transition. Similarly, black chromatophores were retained only if they were observed earlier already as black. Size was determined on hand-drawn ROIs outlining each chromatophore in the aligned dataset.

Colour/Development/numbers model

For the juvenile animal, the birth rate of new chromatophores was estimated by counting, in datasets aligned over days, all chromatophores within a patch of skin on the last day. We then found the fraction of these chromatophores present on previous days within the same aligned patch. The birth rate was calculated as an exponential fit to this data, using 10-fold cross validation of skin regions. The estimated ratio of light/dark was counted from manual annotation of a patches of skin taken at high-resolution (N = 7 animals). In the adult animal, light-to-dark chromatophore transition took longer than our 42-day observation period. The derivation of our model and the method used to estimate light-chromatophore lifetime are provided in full in Supplementary Information.

Growth model

The model described in Fig. 6 is illustrated in more detail in Extended Data Fig. 9a. For simplicity, growth of a skin patch was modelled by a sequence of three steps, repeated every “day”: (i) insertion of new chromatophores, (ii) isotropic expansion of skin patch, (iii) age-dependent size change of inhibitory surrounds and update of chromatophore-colour label. First, we explain how the inhibitory surround was constructed from observations. Second, we define the relevant parameters. Finally, we describe the simulation steps in detail and explain how parameters affect simulation outcome.

The zone of inhibition surrounding individual chromatophores in our spatial model was generated from the empirical average of chromatophore density surrounding a chromatophore (N = 9,199 chromatophores of both colours, 1 animal). We normalized the radial average of the density by the value of the first peak and set any value occurring at greater radial distance to 1. We then fit a logistic function to this density, where x is the distance from the chromatophore center, s the size (i.e. distance at half-height) and k the slope at half-maximum. We inverted this function as 1-I(x), to arrive at the 1D inhibitory surround kernel (Extended Data Fig. 9d). This curve defined the radial dependence of the isotropic 2D inhibitory surround with vector =(x denoting the 2D spatial coordinates.

With the shape of the surround fixed, the simulation contained 5 parameters: 1. The maturation age of chromatophores (L-D transition day); 2. the rate at which chromatophores move away from each other daily (skin growth rate); 3. the size of the inhibitory surround of a chromatophore at birth, s; 4. the rate of change of the inhibitory surround as a chromatophore ages, r; and 5. the threshold level of skin ‘filling’ at which a “day” is complete. Note that these five parameters can all be varied independently of each other. Extended Data Fig. 9a shows the simulation steps at which each parameter is introduced.

Our analytical growth model (above, and Supplementary Information), determines the coupling of maturation age (parameter 1), the rate of new chromatophore insertion (approximately parameter 2 squared), and the L/D ratio. We therefore fixed the values for parameter 1 and 2 approximately to values observed in a young animal and, as expected, observed a realistic L/D ratio. Parameters 3-5 determined the local spatial layout of chromatophores. Parameter 5 allowed us to model the possibility that chromatophores are not inserted to maximum packing density on any given day. Extended Data Fig. 9 illustrates that this parameter affects chromatophore packing without affecting L-D inter-digitation. Our results did not depend on the number of simulated “days”, provided that simulations were run for long enough (~60 “days”) to allow L/D to reach steady state.

We defined an “inhibition”-field F( over the simulated patch of skin of initial size l x l (l = 540 µm), where both components of x= (x are real numbers in the interval The inhibition field was computed as i.e. by adding the inhibitory kernels J (see above) from all present chromatophores i with position and size (determined by parameters 3 and 4). Values larger than 1 were set equal to 1. F was updated after every chromatophore insertion. The simulation began with an empty patch of skin, i.e., with an inhibition-field equal to zero (uniform probability for chromatophore insertion). During a simulation “day” chromatophores were inserted sequentially. A location was drawn from a two-dimensional uniform random distribution over the space covered by the skin patch, and insertion took place with probability p = 1 – F at that location. The “day” ended when no location in the inhibition-field was left with a value less than the filling threshold (parameter 5). At the end of that “day”, the system was expanded by scaling all positions by a fixed rate (parameter 2). The size of the inhibitory surround was then adjusted according to chromatophore age, a, as s, (until a = 19 days, after which s was fixed), and the colour of each chromatophore was updated.

We fit parameters 3-5 to the radially averaged chromatophore-triggered densities from 1 animal (N=4095 D, 5104 L) using a grid search and mean-squared-error loss. Search space was: 3. 67-81 µm radius (full width, half maximum); 4. -3.2%:-1.6%/“day” (adjusted every 2 “days”); 5. 0:0.5 filling threshold.

The parameters of the best fit model were 1: 19 days; 2: 2.06%/“day”; 3: 75.6 µm radius (full width, half max); 4: -2.8%/“day” (adjusted every 2 “days”); 5: 0.1. To approximate squid skin (Extended Data Fig. 10), we adjusted parameter 4 to 2.8%/“day” (adjusted every 2 “days”) until age 45 followed by an expansion to 340 µm radius.

Image manipulation

Colour images in Fig 1c, Fig 1d, Fig 3c, and Fig 5e were uniformly and linearly scaled for clarity.

Statistics

Unless stated otherwise, error bars throughout this study refer to mean ± s.d. For boxplots, central line: median, box limits: quartiles. Whiskers extend to a maximum of ±2.7 standard deviations. No statistical methods were used to predetermine sample size. The experiments were not randomized and the investigators were not blinded to allocation during experiments and outcome assessment.

Supplementary Material

Supplementary Information is linked to the online version of the paper at www.nature.com/nature.