Modeling binary and graded cone cell fate patterning in the mouse retina

Nervous systems are incredibly diverse, with myriad neuronal subtypes defined by gene expression. How binary and graded fate characteristics are patterned across tissues is poorly understood. Expression of opsin photopigments in the cone photoreceptors of the mouse retina provides an excellent model to address this question. Individual cones express S-opsin only, M-opsin only, or both S-opsin and M-opsin. These cell populations are patterned along the dorsal-ventral axis, with greater M-opsin expression in the dorsal region and greater S-opsin expression in the ventral region. Thyroid hormone signaling plays a critical role in activating M-opsin and repressing S-opsin. Here, we developed an image analysis approach to identify individual cone cells and evaluate their opsin expression from immunofluorescence imaging tiles spanning roughly 6 mm along the D-V axis of the mouse retina. From analyzing the opsin expression of ~250,000 cells, we found that cones make a binary decision between S-opsin only and co-expression competent fates. Co-expression competent cells express graded levels of S- and M-opsins, depending nonlinearly on their position in the dorsal-ventral axis. M- and S-opsin expression display differential, inverse patterns. Using these single-cell data, we developed a quantitative, probabilistic model of cone cell decisions in the retinal tissue based on thyroid hormone signaling activity. The model recovers the probability distribution for cone fate patterning in the mouse retina and describes a minimal set of interactions that are necessary to reproduce the observed cell fates. Our study provides a paradigm describing how differential responses to regulatory inputs generate complex patterns of binary and graded cell fates.


Introduction
3 retinoic acid levels, have also been implicated in photoreceptor patterning (Alfano et al., 2011;Satoh et al., 2009). For this study, we focus on modeling the contributions of TH and Thrβ2 to cell fate outcomes.
We desired to quantitatively model cone fate specification in the mouse retina. Our current theoretical understanding of cell fate determination within a tissue describes individual cell types as distinct valleys on an "epigenetic landscape" (Furusawa and Kaneko, 2012;Micheelsen et al., 2010;Wang et al., 2011a;Zhang and Wolynes, 2014;Zhou et al., 2012). Cells make fate decisions by transitioning to one of these "attractor" states on the landscape (Olsson et al., 2016). Differences in gene expression between the states give rise to phenotypic differences between cell types. However, clustering based on single-cell transcriptomics data alone misses subpopulations unless hidden variables are accounted for (Buettner et al., 2015;Setty et al., 2019).
Recently, computational work has also focused on developing mechanistic models of cell-fate decisions (Olariu and Peterson, 2019;Rothenberg, 2019;Teles et al., 2013), especially the formation of patterns in time and space (Formosa-Jordan, 2018;Liang et al., 2015). Multiscale approaches that combine stochastic and deterministic models of tissues at the scale of individual cells have shown promise in helping to elucidate the details of tissue patterning (Coulier and Hellander, 2018;Engblom, 2018;Engblom et al., 2018;Folguera-Blasco et al., 2019;Johnston et al., 2011). The zebrafish retina has been studied to model cell fate decision making based on anticlustering mechanisms that give rise to a lattice structure of differentiated cell types (Cameron and Carney, 2004;Ogawa et al., 2017;Tyler et al., 2005). The near perfect grid of cone subtypes in the zebrafish eye is in stark contrast with the variable patterning of cone subtypes in the D-V axis of the mouse retina (Baden et al., 2013;Haverkamp et al., 2005;Viets et al., 2016).
Here, we present a multiscale model describing the emergence of the complex arrangement of cone cells found in the mouse retina using both stochastic and deterministic methods. We collect data for the model from analysis of immunofluorescence images of retina tissues to identify and map individual cones along the entire D-V axis of the mouse retina. Based on opsin expression in the individual cells, we find that cones can be classified into two main subtypes: S-only cones and co-expression competent (CEC) cones. The S-only cones express S-opsin only, whereas the CEC cones express M-and/or S-opsins in opposing dorsalventral gradients, with higher levels of M-opsin in cones in the dorsal retina and higher levels of S-opsin in cones in the ventral retina. We then use the data to parameterize a mathematical model of a two-step cone patterning process.
Step one is a binary choice between S-only fate and CEC fate. If CEC fate is selected, a second mechanism regulates S-and M-opsin expression in a reciprocal, graded manner, along the dorsal-ventral axis. Our quantitative modeling shows that the expression of S-and M-opsins in CEC cells are differentially activated based on dorsal-ventral patterning inputs from T3. Our model closely recapitulates cone patterning observed in the mouse retina and provides insights into how spatial patterning inputs regulate binary and graded features of cell fate in parallel.

Characterization of cone subtype patterning in the mouse retina
To globally characterize patterning of opsin expression in the mouse retina, we first examined the relative intensity of S-and M-opsin expression in the D-V and temporal-nasal (T-N) axes at low resolution for whole-mounted retinas. We immunostained and imaged S-and Mopsin proteins at 100X magnification (see Materials and Methods). Following image acquisition, we manually rotated each image so that the D-V axis was aligned vertically (Fig. 1A, E, I). At this resolution, individual cells cannot be identified, so we instead subdivided each image using a 25 pixel x 25 pixel grid, which is an area containing approximately one to two cells. Within each bin of the grid, we counted the number of pixels that had significant S-opsin signal alone, M-opsin signal alone, or both M-and S-opsin signals. We then normalized each bin by the total number of pixels with expression in that bin. This calculation gave us the relative density of each photoreceptor type by location in the retina (Fig. 1B, F, J).
Next, we quantified global differences in patterning in the D-V and temporal-to-nasal (T-N) dimensions. We averaged the binned density values to obtain the relative density as a function of either D-V (Fig. 1C, G, K) or T-N position (Fig. 1D, H, L). We observed distinct transitions in both S-and M-opsin expression along the D-V axis. High levels of M-opsin in the dorsal region exhibit a gradual transition to low levels in the ventral region (Fig. 1C). In contrast, S-opsin shows a rapid transition from zero to high expression in the D-V axis (Fig. 1G). As these opsins display an inverse yet non-complementary relationship, co-expression was most prominent in the middle third of the retina where these two transitions overlap (Fig. 1K). We observed minimal variation in S-and M-opsin signal in the T-N axis (Fig. 1D, H, L). We imaged and analyzed six wild-type retinas at this resolution and saw a similar pattern in each (Table  S1). Together, we observed differential graded patterning for S-and M-opsin expression along the D-V axis (Fig 1C, G, K).

Analysis at single cell resolution reveals two distinct cone subtype populations
To further investigate photoreceptor patterning along the D-V axis, we next analyzed cone subtype specification at the single-cell level. For the same six retinas, we imaged S-and M-opsin expression at 200X magnification in a strip measuring approximately 600 m x 6,000 m aligned vertically along the D-V axis ( Fig. 2A). Previous studies analyzed ~500 m in the dorsal-ventral axis centered on the transition region (Haverkamp et al., 2005), whereas our approach enabled evaluation of the entire ~6,000 m length of the retina.
At 200x magnification, we were able to distinguish and identify individual cells. We developed an analysis pipeline to identify the position, size, and boundaries of each cell, a process known as segmentation (see SI Methods, 1.1.1). Overall, we identified ~250,000 total cells across six retinas. Using these cell boundaries, we calculated the expression intensity of M-and S-opsin for each cell. We classified cones into groups expressing M-opsin only (Fig. 2E, K, Q), S-opsin only (Fig. 2F, L, R), and S-and M-opsin co-expression (Fig. 2G, M, S). The pipeline did not identify distinct morphological or size differences among the cones (data not shown). We found that the pipeline's accuracy and false positive rate were comparable to handscored retinas (see SI Methods 1.

1.2).
After obtaining the cell boundaries, we quantified the density of cone subtypes based on opsin expression relative to D-V position. Consistent with our low-resolution analysis, we observed a gradual decrease in the abundance of cones expressing M-opsin in the dorsal to ventral direction (Fig. 3A), contrasted by a sharp increase in S-opsin expressing cones (Fig.  3B). We fit these curves to Hill functions to quantify the steepness of the transition (Fig. S1). The transition in the S-opsin expressing cells is extremely sharp with an average Hill coefficient of ~30 while the M-opsin transition is much more gradual with a coefficient of ~2-3.
To compare the transition region between retinas, we established a reference point to align the images. Since the S-opsin transition is sharp and an external reference is absent, we used the midpoint of the S-opsin transition from the fit as the reference point. We aligned all of the retinas and overlaid the transition fits (Fig. S2). The relative position of the S-opsin and Mopsin transitions are consistent from retina to retina, suggesting that the transitions in S-opsin and M-opsin expression are driven by a common effector. M-opsin only expressing cones decline at the transition point (Fig. 3C), coincident with the dramatic increase in S-opsin expression (Fig. 3B). At this transition point, cones begin to express both S-and M-opsins (Fig.  3E) and the cone populations are very diverse, comprised of those expressing M-only, S-only, and varying levels of both S-and M-opsins ( Fig. 2H-M). The fraction of S-only cells gradually increases from ~1% of cones in the dorsal region to ~20-30% in the ventral region (Fig. 3D, S3).
These analyses show the differential, inverse responses of S-and M-opsin expression to D-V patterning inputs on the individual cell level.
In a previous study, Haverkamp et al. measured differential opsin expression of cone cells in a window of ~500 m near the transition point (Haverkamp et al., 2005). In agreement with our data, they observed that ~8-20% of cones expressed only S-opsin. Our results show that this measurement was part of a broader binary decision trend extending much further along the D-V axis in both directions. They also discovered that within this population, in the ventral region where S-only cones are more abundant, about 5% of S-opsin only cones contact S-cone bipolar cells and they classified these as genuine S-cones. These genuine S-cones are evenly distributed across the retina (Haverkamp et al., 2005).
To see if we could distinguish these two classes of cone subtypes in our data, we calculated the joint probability distribution for S-and M-opsin intensity in individual cells for each retina (Fig. 4A, S4). When considering all cones in the retina, there appear to be three clusters of cell-types corresponding to the three classifications that we defined earlier: S-only, M-only, and co-expressing. However, when the data are analyzed by D-V position ( Fig. 4B-F, S5), a different pattern emerges. First, we see an S-only cluster than has high and consistent expression of S-opsin while increasing in abundance along the D-V axis (Fig. 4). Interestingly, the other two clusters resolve into a single cluster that changes position in a continuous way, gradually moving from low S-opsin expression and high M-opsin expression in the dorsal region, to moderate S-opsin expression and low M-opsin expression in the ventral region (Fig. 4).
Thus, these data suggest that the mouse retina contains two main subtypes of cones: 1. S-only cones that have high S-opsin expression independent of D-V position and 2. coexpression competent (CEC) cones that express S-and/or M-opsins dependent upon D-V position. In the ventral region there is a mixture of S-only cones and CEC cones that express Mopsin at a very low level. It is difficult to distinguish these two classes using only S-opsin expression, but as can be seen in Fig 4., the two populations are well separated when comparing both M-and S-opsin intensities. The S-only cones identified with our approach may, therefore, contain a subclass corresponding to the genuine S-cones identified by Haverkamp et al., but as we could not distinguish their connectivity to bipolar cells, we are only able to describe the populations of S-opsin only expressing cones.

Expression levels of S-and M-opsin in cone cell subtypes
Having classified the major subtypes of cones and related their positions and opsin expression states, we next evaluated the D-V dependence of the opsin expression intensity in individual cones. We quantified opsin expression for all M-opsin expressing cones (Fig. 5A, F), all S-opsin expressing cones (Fig. 5B, G), M-and S-opsin CEC cones (Fig. 5C, H), M-opsin only CEC cones (Fig. 5D, I), and S-opsin only cones (Fig. 5E, J) relative to their D-V positions within the retina.
In CEC cones, M-opsin expression levels decrease in the D-V axis, with the midpoint of expression level located at the transition point ( Fig. 5A, C, D, S6). In contrast, S-opsin expression in CEC cells is very low in the dorsal region and increases linearly in the D-V axis starting at the transition point ( Fig. 5F, H, I, S6). The slope of increase for S-opsin is steeper than for the M-opsin decrease (Fig. S7).
Compared to CEC cones, S-only cones have an overall higher expression level of Sopsin, particularly in the dorsal region (Fig. 5J compared to G and H). M-opsin expression in Sonly cones is significantly lower than the lowest M-opsin expression seen in CEC cones (Fig.  5E). In Fig. 5B, this difference can be seen as two distinct lines of density (Fig. 5B, arrow heads). Together, these analyses defined the expression of S-and M-opsin in the two cone populations in relation to their D-V positions in the retina.

Modeling cone subtype fate decisions
To interrogate how regulatory inputs could produce the complex pattern of binary and graded cell fates in the mouse retina, we developed a multiscale model describing the probability distributions of the cone subtype decisions (i.e. binary choice) and S-and M-opsin expression levels (i.e. graded) as functions of position along the D-V axis (Fig. 6, SI Methods 1.2). We modeled a 5 mm x 1 mm x 5 µm section of the retina with the long dimension aligned with the D-V axis.
TH signaling activates M-opsin expression and represses S-opsin expression (Ng et al., 2001;Roberts et al., 2006). T3 is a critical regulator of cone subtype fate in the human retina (Eldred et al., 2018), and scRNA-seq data suggest that Thrβ2 is expressed in all mouse cones (Clark et al., 2019). Though other diffusible factors and transcription factors play roles (Alfano et al., 2011;Roberts et al., 2005;Satoh et al., 2009), TH signaling is the main and best-understood determinant of cone subtype fate. Thus, we built a simplified model of cone subtype specification based on the dorsal-ventral regulation of cone fates by the gradient of T3.
Within the modeled volume, T3 molecules diffuse according to the deterministic diffusion equation with constant concentration boundaries, establishing a D-V gradient (Eq. S1, SI Methods). Also, within the volume, we stochastically modeled ~23,000 individual cones spaced on a hexagonal grid (Fig. 7). These cells stochastically exchange T3 molecules with the surrounding deterministic microenvironment. Within each cone, T3 can bind to and activate Thrβ2 (Thrβ2*), controlling both fate specification and opsin expression ( Fig. 6B1-4).
For the binary fate decision, we defined a fate determinant function, FD(X). Photoreceptors start in an undifferentiated fate, FD(U), and stochastically progress to either the FD(S) (S-only) fate or FD(C) (CEC) fate (Fig. 6B2). Selection of the FD(S) fate is negatively influenced by Thrβ2* (Fig. 6B2). Once cells enter the FD(S) fate, S-opsin is constitutively expressed at a high level regardless of D-V position (Fig. 6B2). In FD(C) cones, M-opsin expression is induced by Thrβ2* (Fig. 6B3). Conversely, S-opsin in FD(C) cells is negatively regulated by Thrβ2* and positively regulated by inactive Thrβ2 receptors (Fig. 6B4). Full details of the model are given in the SI Methods 1.2, along with parameterization details.

Model-based simulations recapitulate experimental cone patterning
To compare the output of our probabilistic model to experimental data sets, we ran a set of 100 individual simulations and calculated the probability distributions of various observables. Fig. 7 shows the output of one simulation. Moving from dorsal to ventral, the model reproduces the gradual increase in the fraction of S-only cones, FD(S) (Fig. 7A1-4), as well as the sharp transition in CEC cones, FD(C), expressing S-opsin at the transition zone ( Fig. 7B1-4). Similarly, we observed the gradual decrease in the fraction of CEC cones expressing M-opsin ( Fig. 7C1-4). In the overlapping region, there are a significant number of cones that co-express both S-and M-opsins ( Fig. 7D1-4). In the dorsal region, a small number of S-only cones that highly express S-opsin are readily apparent ( Fig. 7A1-2, E1-2).
To characterize how well our model recapitulated the observed experimental cell distributions, we calculated the mean density of cells of various phenotypes as a function of D-V position (Fig. S8). These average density profiles compare well to the experimental density profiles shown in Fig. 3. Together, cone fate patterning and expression levels are highly similar in our model and the imaged retinas: S-only cells ( We parameterized our model using the mean of all the retinas sampled, which exhibited retina-to-retina variability (Fig. S1, S3, S7). Therefore, it is not expected that our model will exactly recapitulate the patterning of any individual retina.
We next calculated the probability distributions of S-and M-opsin expression along the D-V axis for our simulation data (Fig. S9). The mean intensity of M-opsin in CEC cones gradually decreases as D-V position increases. The S-opsin distribution shows high expression in the ventral-most region, but has two separate populations in the dorsal region: the highly expressing S-only cells and the lowly expressing CEC cells. The CEC cones converge to zero S-opsin expression in the dorsal region while the S-only cones maintain high expression as they decrease in abundance. Because our simulated distributions are constructed from 100 independent simulations, the probability density of the S-only cones is much smoother than in the experimental data (compare Fig. S6 and S9). The simulated expression features are in agreement with the experimental expression profile (Fig. 5, S6).
We next related the joint probability distributions for the experimental (Fig. 8A-E) and simulated ( Fig. 8F-J) data along the D-V axis. The simulated and experimental data show two distinct populations: 1) S-only cones with high S-opsin expression and no M-opsin expression whose expression levels are independent of D-V position, and 2) CEC cones that gradually change from high M-opsin and low S-opsin expression to moderate M-opsin and high S-opsin expression along the D-V axis. Fig. S10 shows the joint probability distribution between S-and M-opsin expression within 250 µm D-V bins for all 100 simulations. In the high resolution simulated data, it is evident that the position of the CEC cell cluster gradually changes with D-V position. Our model closely simulates the experimental data, and supports the hypothesis that cells respond differentially to the same morphogen gradient, producing both binary and graded cell fates.

Correlation between S-opsin and CEC fate decisions
Our quantitative simulations give us the capacity to test various hypotheses about retinal patterning. We wanted to know whether the gradual decrease in the CEC cone population and the sharp increase in S-opsin expression in these CEC cells were driven by a shared upstream signaling input. If these two processes respond to the same upstream input, we would expect that they should be coupled and be linked by D-V position. If, however, they do not respond to the same input we would expect that their transitions should be independent. In our model, they are coupled through the T3 gradient and we wanted to test if the experimental retinas displayed more or less variability than expected. Because stochasticity in both the experimental data and our model introduces retina-to-retina variability, we checked for overlap of the corresponding probability distributions.
We calculated the probability of having a given CEC population fraction at the S-opsin transition point from both our simulations and experimental retinas (Fig. 9A). We also calculated the rate at which the CEC population decreases at the transition point (Fig. 9B). As can be seen from both plots, the two probability distributions are in good agreement. In particular, the widths of the experimental probability distributions are similar to the widths from the stochastic simulations. With the small number of experimental data points, we do not assign a level of statistical significance to the overlaps, but they provide qualitative evidence that the two transitions are in fact coupled through a shared upstream input.

Comparison of Thrβ2 mutant retinas to wild-type retinas
Finally, to elucidate the effect of the T3 gradient on cone cell patterning in the mouse retina, we dissected and imaged Thrβ2 knockout mutant retinas (ΔThrβ2). In the absence of functional Thrβ2, no M-opsin is expressed (Applebury et al., 2000;Ng et al., 2001;Roberts et al., 2006). Consistent with previous work, we observed no M-opsin expression in these mutant retinas. Fluorescence from anti-M-opsin antibodies was nonspecific and stained cell and background with equal intensity (Fig. S11).
The expression of S-opsin in cones was also markedly different between ΔThrβ2 and WT retinas. Both the density of S-opsin-expressing cones and the expression distribution are flat with respect to the D-V axis (Fig. 10B, S12). Also, the relative intensity of the S-opsin signal across the retina was much lower in ΔThrβ2 retinas than the maximum value seen in WT retinas (e.g., from cones in the ventral region). ΔThrβ2 retinas and wildtype retinas were taken at the same time and stained with the same batch of antibody, then imaged with the same laser intensity for comparison of opsin levels. We found that the relative intensity of opsin staining in cones was most similar to the middle region of WT retinas. This effect is consistent with our model in which S-opsin expression is controlled through a combination of negative regulation by active Thrβ2* and positive regulation by inactive Thrβ2 (Fig. 6B4).

Discussion
In these studies, we described the distribution of cone photoreceptors in the mouse retina and developed a quantitative model for the specification of binary and graded cell fates in response to D-V regulatory inputs. By using high-resolution microscopy combined with automated image analysis, we expanded on previous studies and mapped the cell fate decisions of cone cells across an entire dorsal to ventral region of the mouse retina. By analyzing cell fates in the context of their position in the tissue, we found that cones could be classified into two subclasses with a graded gene expression profile changing nonlinearly. This study exemplifies the benefits of quantitatively analyzing populations of cells in a tissue when classifying fate decisions.
In the mouse retina, we defined two cone subtypes, S-only cones and CEC cones, based strictly on opsin expression profiles. Interestingly, the population of S-only cones in the dorsal region have higher S-opsin expression than most S-only cones in the ventral region (Fig.  5J). These highly expressing S-only cones are found at a steady density across the D-V axis of the retina. It is possible that this subset of evenly distributed, high S-opsin expressing cells could comprise the "genuine" S-cones that connect to blue-cone bipolar neurons (Haverkamp et al., 2005). We were not able to mathematically distinguish genuine S-cones from the total population of S-only cones. Together, opsin expression and connectivity suggest three possible distinct cone subtypes: 1. CEC cones that do not connect to blue-cone bipolars, 2. "genuine" Sonly cones that connect to blue-cone bipolars, and 3. S-only cones that do not connect to bluecone bipolars. Examination of connectivity in conditions that perturb thyroid hormone signaling may inform the relationship between opsin expression and connectivity and their relationship to cone subtype.
We developed a mathematical model that described both the binary fate specification process of cones and the graded expression of opsins, all driven by an external gradient. Stochastic modeling of this complex process generated probability distributions that we used to compare with experimentally observed cell distributions to test hypotheses about the connections between cell fates. Stochastic modeling is now sufficiently mature to perform detailed simulations of tissue-level cell-fate decisions. Combining stochastic models with highthroughput microscopy is a powerful tool for helping to understand complex relationships in tissues.
These methods advance our understanding of how regulatory inputs influence complex cellular decisions to specify binary and graded cell fates within the same cell type in the same tissue. A next step will be to integrate more signaling inputs into our model for retinal development. Numerous signaling molecules are expressed in D-V gradients and are involved in retinal and cone cell development. Specifically, retinoic acid (RA) is an important morphogen that is expressed at high levels ventrally during development, and then at moderate levels in the dorsal region in the adult mouse (McCaffery et al., 1992;McCaffrery et al., 1993). Moreover, further studies would include integrating transcription factor binding partners of Thrβ2 into the model, as Thrβ2 acts as a homodimer and as a heterodimer with RXRγ (Roberts et al., 2006). This work represents an important first step towards modeling the complex network of interactions that guide binary and graded cell fate specification.
The retina provides an excellent paradigm to study how signaling inputs generate patterns in two dimensions. The next challenge will be developing models for patterning in more complex 3-dimensional neural tissue found in brain structures. Quantitative modeling has enormous potential to integrate multiple signals across a tissue and build networks to better understand and predict the outcomes of development when variables are changed, for instance in disease states.

Materials and Methods: Animals
Mice (strain C57BL/6) were housed under a 12 h light:12 h dark (T24) cycle at a temperature of 22°C with food and water ad libitum. Male and female mice at 2-8 months of age were housed in plastic translucent cages with steel-lined lids in an open room. Ambient room temperature and humidity were monitored daily and tightly controlled. Wild-type mice (C57BL/6; Jackson Laboratory), and Thrb tm2Df mutant mice (gift from the Forrest Lab) were used in this study. Thrb tm2Df mutant mice specifically knock out expression of Thrβ2, and leave Thrβ1 intact as previously described (Ng et al., 2001). All animals were handled in accordance with guidelines of the Animal Care and Use Committees of Johns Hopkins University. All efforts were made to minimize the pain and the number of animals used.

Immunohistochemistry
Retinas were dissected in PBS, then fixed in fresh 4% formaldehyde and 5% sucrose in PBS for 1 hour. The dorsal portion of the retina was marked with a cut. Tissue was rinsed 3X for 15 min in PBS. Retinas were incubated for 2 hours in blocking solution (0.2-0.3% Trition X-100, 2-4% donkey serum in PBS). Retinas were incubated with primary antibodies in blocking solution overnight at 4°C. Retinas were washed 3X for 30 min in PBS, and then incubated with secondary antibodies in blocking solution for 2 hours at room temperature. At the end of staining, retinas were cut to lay flat on a slide, and were mounted for imaging in slow fade (S36940, Thermo Fisher Scientific).

Microscopy and image processing
Fluorescent images were acquired with a Zeiss LSM780 or LSM800 laser scanning confocal microscope. Confocal microscopy was performed with similar settings for laser power, photomultiplier gain and offset, and pinhole diameter. Whole retinas were imaged with a 10X objective, and maximum intensity projections of z-stacks (5-80 optical sections, 4.9 μm step size) were rendered to display all cones imaged in a single retina. Retinal strips were imaged with a 20X objective, and maximum intensity projections of z-stacks (5-80 optical sections, 1.10 μm step size) were rendered to display all cones imaged in a single retina. ΔThrβ2 retinas and wildtype retinas were taken at the same time and stained with the same batch of antibody, then imaged with the same laser intensity for comparison of opsin levels.

Segmentation of cone cells from microscopy images
Microscopy images were analyzed using a custom parallel image processing pipeline in Biospark (Klein et al., 2017). Briefly, each fluorescence channel was first normalized and filtered to remove small bright features. Then, each remaining peak in fluorescence intensity was identified and an independent active contour segmentation (Marquez-Neila et al., 2014) was performed starting from the peak. If the resulting contour passed validation checks it was included in the list of segmented cone cells for the channel. Finally, cell boundaries were reconciled across both channels to obtain a complete list of identified cone cells. Full details are given in the SI Methods.

Modeling cone cells fate decisions and opsin expression
Multiscale modeling of the retina strip was performed using a hybrid deterministic-stochastic method. Diffusion of T3 in the microenvironment of the retinal strip was modeled using the diffusion partial differential equation (PDE). The PDE was solved using an explicit finite difference method. The cone cells were modeled using the chemical master equation (CME) to describe the stochastic reaction scheme implementing the cell fate decision-making. The CME for each of the 23,760 cone cells was independently sampled using Gillespie's stochastic simulation algorithm (Gillespie, 1977). Reconciliation between the CME trajectories and the PDE microenvironment was done using a time-stepping approach. Complete mathematical details of the model and simulation methods are available in the SI Methods. All simulations were performed using a custom solver added to the LMES software (Roberts et al., 2013), which is available on our website: https://www.robertslabjhu.info/home/software/lmes/.

Segmentation of individual photoreceptor cells in 200X images
To identify and segment individual cone cells in the immunostained fluorescence images, we developed an image processing pipeline similar to that used previously [1]. The analysis begins Next, we find the outline of each potential cell using an active contouring method known as morphological snakes [2]. Because we are segmenting tens-of-thousands of individual potential cells from each image, we perform this step of the pipeline in parallel using the Biospark framework [3], which is a data intensive parallel analysis package for Python. We perform validation of the segmented boundaries before classifying the object as a cone cell. https://www.robertslabjhu.info/home/software/mouse_eye.

Validation of segmentation results
We validated that our segmentation algorithm produced results similar to human annotators by comparing manually and automatically generated statistics from representative samples of our data set. We picked seven different regions and manually counted the number of S-opsin only, Mopsin only, and coexpressing cells. We then analyzed the same regions using our segmentation algorithm and obtained the automatically generated classifications. Table S2 shows the counts of cell types from these regions from both human and computer annotations. Overall, there is excellent agreement in the relative abundance of the different cell types. The absolute counts have some systematic difference, with the automatic segmentation typically identifying more cells than human annotators. This is mostly due to what appear to be single long cells that are split into multiple bright pieces separated by low fluorescence breaks.
Human annotators tend to regard the trace as a single long cell, while the automatic segmentation tends to identify multiple smaller cells. Importantly, we do not know the true underlying cell morphology, so we cannot generally say whether the human or automatic annotators are more accurate. In any case, since the cell parts are identified correctly and the segmentation is consistent across retinas, we expect these minor difference to have no impact on our results.

Modeling cone cell fate determination and opsin expression in a retinal strip
To model the cone fate decisions in a large retinal strip we use a combination of stochastic and deterministic modeling. We start with a three-dimensional volume 5 mm long in the X dimension, 1 mm wide in the Y dimension, and 5 µm in the Z dimension representing the microenvironment of the dorsal-ventral (DV) strip. The small z dimension make this an effectively two-dimensional system and we include z below only for completeness. Within this volume we model diffusion of thyroid hormone (T3) using the deterministic diffusion equation: where C(r, t) is the concentration of T3 at position r and time t, D is the diffusion coefficient used for T3, and ∇ 2 is the Laplace operator ( ∂ 2 C ∂x 2 + ∂ 2 C ∂y 2 + ∂ 2 C ∂z 2 ). See Table S3 for all parameter values used.
We numerically solve the diffusion partial differential equation (PDE) using a explicit finite difference method with grid spacing dx and a time step dt = (dx 2 )/(2 · 6D), where the extra factor of 2 in the denominator ensures numerical stability. We fix the concentration at the X = 0 boundary to C hi and at the opposite boundary to C lo to establish a stationary concentration gradient in the X dimension. The Y and Z boundaries are taken to be reflective. We initialize the concentrations C(r, 0) according to a linear decrease in r to follow to the boundary conditions. Within the microenvironment, we place 23,760 individual photoreceptor cells spaced on a hexagonal grid spanning the X-Y plane, with center-to-center distance d cell−cell and with a radius r cell . Each cell is modeled independently using the chemical master equation (CME) describing the probability to have a particular count of each species: where P t (x) is the probability for a cell to have a particular state vector x giving the count for each chemical species and A is a transition matrix describing all of the reactions between the chemical species.
Within each photoreceptor cell a series of reactions describing the fate of the cell and also opsin expression take place. First, cells contain thyroid hormone receptors THRβ2. T3 can bind reversibly to THRβ2 to switch it to an activated state THRβ2 * : Note that T3 also diffusing across the PDE microenvironment and the value of T3 is synchronized between the PDE and CME models. The synchronization procedure is discussed below.
Next, the cells switch between three cell fates. Little is known about these steps, but the fate decisions appear to be stable. Therefore, we model photoreceptor fate decision-making as barrier crossing process with a variable number of cooperative steps n.
The FD(U) ↔ FD(S) transition is described by Here, the syntax FD(U→ i S) denotes a cell that has progressed i steps along the path from the FD(U) to FD(S) fate. Also, H is an inhibiting Hill-like kinetic function defined by with k lo and k hi the lower and upper limits of the kinetic process, respectively, k m the midpoint of the transition, and h the Hill exponent giving the cooperativity of the transition. Likewise, the FD(U) ↔ FD(C) transition is described by Then, both S-opsin and M-opsin proteins can be expressed by photoreceptor cells, depending on the cell type and the local concentration of T3. We model opsin expression using the following kinetic equations H is an activating Hill-like kinetic function (S14) Finally, both opsins can be degraded according to We initialize each cell at t = 0 to the FD(U) state with a copy number of THRβ2 proteins independently sampled from a Gaussian distribution with a mean concentration of µ thrb and a variance of σ 2 thrb . The number of T3 molecules available to the cell is initialized according to the T3 concentration at the cell's D-V position. Likewise, the fraction of activated THRβ2 is initialized to its equilibrium value according to the cell's D-V position. All opsin counts are initialized to zero.
We then model the stochastic time evolution of each cell using the standard Gillespie stochastic simulation algorithm (SSA) [4,5].

Microenvironment modeling of combined PDE and CME dynamics
Since we are performing two parallel simulations, PDE and CME, we need to partition the molecules between them. Each cell has a volume smaller than a PDE subvolume and each cell is assumed to be completely contained within a single subvolume. We initialize the T3 molecule count in the extracellular space of each cell to be the rounded number of molecules corresponding to the cell's extracellular volume multiplied by the PDE subvolume's concentration. In this way T3 molecules are represented in each simulation.
To integrate the PDE and CME dynamics, we implemented a parallel time-stepping approach.
We divide time into discrete synchronization intervals ∆t and evolve overall time according to t i+1 = t i + ∆t. (S17) During each ∆t we update the state of each cell using the SSA and the concentration of each PDE subvolume using a finite difference algorithm. The loss or gain in the number of extracellular T3 molecules in the stochastic simulation during the time step is tracked. This gain or loss is converted into a concentration flux and applied to the PDE subvolume during the next time step.
Simultaneously, the number of extracellular T3 molecules in the stochastic model is reset for the next time step according to the new subvolume concentration. Our method conserves mass and has good error characteristics as long as the T3 flux at each time step is of the order of a few molecules. The method is implemented as the "microenvironment" solver in our LMES software and is freely available on our website: https://www.robertslabjhu.info/home/software/lmes.  (5) the per cell S-opsin expression level. Because our data were collected from multiple retinas, we first fit the raw data to functions that we could use to describe a hypothetical mean retina.

Parameterization of retinal strip microenvironment model
We describe the fraction of cells in the various subfates as (modified) Hill-like functions. The fraction of cells expressing S-opsin as a function of D-V position x is described by: where m and b are the slope and x-intercept of a baseline fraction, respectively, x mid is the midpoint of the transition, and h is the Hill coefficient. Likewise, the fraction of M-opsin expressing cells is given by: and the fraction of FD(S) cells is given by: Figures S1+S3 show the fits of these functions to the raw data for the various retinas. We then took the mean of the various parameters to construct a hypothetical mean retina. Figure S13 shows the fraction of cells in these states as a function of D-V position in our mean retina.
We describe the mean per cell expression level of M-and S-opsin as piecewise linear functions with a low and a high limit separated by a biphasic region with two different slopes: Here, y 0 and y 1 are the left and right baselines, x 0 , x 1 , and x 2 are the D-V points where the slope changes, and m 1 and m 2 are the two slopes. Figure S7 shows the fit to the M-and S-opsin expression for the various experimental retinas. Figure S13 shows the values of these functions for our hypothetical mean retina.
Once we had a hypothetical mean retina, we used it to parameterize our model. We derived expressions for the mean value of the various observables as a function of D-V position by solving the deterministic system of ordinary differential equations described by Equations S3-S16. We then used non-linear least squares with Nelder-Mead minimization to globally optimize the parameters using all five data sets. Table S3 gives the best fits values for the free parameters and Figure S14 shows a comparison of the deterministic model calculated using the best-fit parameters against our hypothetical mean retina. Retina Filename Genotype Num. Cells