Procedural generation of synthetic retinal vasculatureAlgorithms for generating synthetic retinas followed multiple, length scale dependent steps66. Firstly, the values of geometrical parameters were set by sampling from a normal or uniform distribution according to parameter values shown in Supplementary Table 1. Retinal vasculature took the form of spatial graphs (i.e. branching nodes connected via vessel segments) which were initially seeded using a Lindenmayer system (L-system), and then extended using constrained constructive optimisation (CCO) and Lattice Sequence Vascularisation (LSV). For the L-system (see Github directory: retinasim/lsystemlib/) retinal artery and retinal vein segments were initially directed from a central optic disc point and branched by five generation to form seeding retinal arcades. These were passed into the CCO algorithm, followed by LSV, to further branch the vessels in progressively finer resolutions. Once complete, the macula region was cleared and regrown into a radially-spaced lattice (https://github.com/AndrewAGuy/vascular-networks). Finally, spatial undulations were imposed onto the vessels by overlaying sinusoidal curves with randomised frequencies and amplitudes (retinasim/sine_interpolate.py). These implementation of these linked algorithms is described further in the sections below.Lindenmeyer system seedingFirstly, seeding networks following approximate retinal vascular branching geometry were constructed, starting with a putative central retinal artery and retinal vein positioned at the centre of the optic disc. The diameters of the retinal artery and vein were 135 ± 15 μm and 151 ± 15 μm, respectively35, oriented parallel to the optic nerve (defined as the z-direction). Two branches were added to the end of each of these segments, oriented in the x-y plane, and one directed above and the other below the retinal midline. Subsequent branching of these vessels was performed stochastically, with segment lengths between bifurcations set as a fixed fraction of vessel diameter (18 ± 3) and bifurcation vessel diameters set according to:$$\cos {\theta }_{1}=\frac{{\left(1+{\alpha }^{3}\right)}^{\frac{4}{3}}+{\alpha }^{4}-1}{2{\alpha }^{2}{\left(1+{\alpha }^{3}\right)}^{\frac{2}{3}}}$$
(1)
$$\cos {\theta }_{2}=\frac{{(1+{\alpha }^{3})}^{4/3}+1-{\alpha }^{4}}{2{(1+{\alpha }^{3})}^{2/3}}$$
(2)
Normally-distributed noise was added to branching angle values, with a standard deviation of 5°. Vessel bifurcation angles were assigned such that the larger vessel oriented towards the macula to create putative major vessels oriented around the macula. Fifth-order bifurcations were added to the network, or until vessels breached the edge of the retina domain. Supplementary Fig. 1a shows an example of an L-system seeding network.Major vessel growthSeed vessel networks were used as input into a multi-scale growth algorithm for the creation of hierarchical vasculature. First, seed networks were amended to provide a uniform distribution of leaf nodes (terminating arteriole and venuole nodes created prior to construction of capillary networks at a later stage) throughout the circular domain, using Accelerated Constrained Constructive Optimisation40, using a leaf node spacing of 3 mm. Multiscale, two-dimensional lattices were defined (stride lengths ranging from 3000 to 150 µm, with five iterations linearly spaced within that range) and used to grow vessel networks by progressively adding vessels into unoccupied lattice sites from neighbouring occupied sites, choosing the candidate vessel which minimised the expected change in network cost (see below), and progressively reducing the length scale when no more progress could be made. After the initial growth stage, all existing leaf nodes were removed67. At all stages of the major vessel growth the macula region was kept free of vessels by removing vessels which intersected it, forcing flow to divert around it.As retinal vasculature is positioned in front of the retina itself, we optimised networks to minimise the area of the retina occluded by vessels, according to a cost function based on Murray’s law:$$C\left(B,\,\lambda,\,\rho \right)={\sum}_{b\in B}\,{r}_{b}^{\rho }{l}_{b}^{\lambda }$$
(3)
with ρ = 1 and λ = 1, and where B is the set of vessel segments in the network, with length l and radius r. After each growth step the network geometry was optimised by moving vessel nodes, and highly asymmetric bifurcations were trimmed for regrowth40 using the thresholds from ref. 39 to account for the high asymmetry of optimal networks39. After growth at each length scale was terminated, the networks were optimised topologically by allowing asymmetric bifurcations to move their low-flow side downstream and branches which were short compared to their expected length under the West, Brown and Enquist model68 to be treated as a single higher-order split for regrouping using a method similar to69. Due to the two-dimensional nature of the networks, network self-intersections were tested for using the approach of40 however, rather than resolving the intersections by making excursions around the contact site we rewire the vessels to prevent future iterations from recreating the same intersection.Unlike the implementation of ref. 40, leaf nodes were allowed to move from their nominal location up to a specified “pinning distance”, given as a fraction of the leaf spacing. Existing vessels could be specified as frozen, in which case the optimiser did not touch them. This approach was used to perturb the optimal root vessel structure with artificial tortuosity, strip away the downstream branches and regrow the downstream vessels, repeating this down the tree structure.Macula growthVessels supplying the macula have a characteristic radial structure, motivating the development of a particular approach to enforce this structure. This uses the same lattice site invasion approach between the macula outer radius and the fovea (which is kept vessel-free), but with the stride set low enough that the majority of the growth arises from spreading over many iterations at the same length scale rather than hierarchical refinement. The macula has a configurable flow rate density compared to the rest of the retina, ranging from 1.5 – 2.0 and leaf nodes are offset by uniformly sampling an offset in a disc around the nominal position to ensure that vessels did not align along the lattice sites. The macular vessels were prevented from doubling back on themselves by setting a hard limit on the vessel angle, preventing obviously non-physiological structures from arising whilst still allowing the radial pattern to develop. After all leaf nodes are created, a sparsity factor is specified and each leaf node removed with this probability, then the remaining vessels are geometrically optimised.Network overpass and interleavingIn the final stage, the arterial and venous networks have their collisions resolved using the method of ref. 40, creating out-of-plane excursions around contact sites between the networks. To enable further micro-scale network growth techniques to create an interdigitated structure, we remove the low-flow side of all arterio-venous intersections with a radius below a critical value (5 um), leaving surviving vessel geometry untouched. Interdigitations were then created using a Space Colonisation implementation70, interspersed with geometric optimisation.Vessel tortuosityThe multi-scale growth algorithm creates relatively straight paths between branching points, and to simulate tortuous retinal vessels, particularly in veins, sinusoidal displacements were overlaid. Two oscillations were superimposed according to:$${d}^{{\prime} }\left(x,\,r\right)=d\left(x,\,r\right)+{a}_{0}\sin \left(\frac{x}{{\tau }_{0}\left(r\right)}+{\delta }_{0}\right)+{a}_{1}\sin \left(\frac{x}{{\tau }_{1}\left(r\right)}+{\delta }_{1}\right)$$
(4)
where d(x,r) is the path taken by a vessel with radius r, and d’(x,r) is the modulated path. The amplitude of displacements, a0 + a1 ranged from r –3.5r for arteries and r –7.5r for veins, with a low frequency period (τ0, ranging from 15r – 25r) and a high frequency period (τ1, ranging from 30r – 50r). The phase of the modulations, φ0 and φ1, enabled modulations to be matched between vessel bifurcations.Simulating vascular flow and fluorescein deliveryBlood flow in retinal networks were simulated using our REANIMATE platform47, which uses a connectivity-based formalism to optimise Poiseuille flow in tree-like spatial graphs. As anterior retinal vasculature features a single arterial inlet and venous outlet, the system requires only one pressure boundary condition (the difference between arterial and venous inlet pressures), which was fixed at 56.2 ± 14.0 and 20.0 ± 10.0 mmHg, respectively.Time-dependent delivery of contrast agent (e.g. fluorescein) was simulated as described in d’Esposito et al.47. Briefly, a bolus of fluorescein was simulated according to$$C\left(t\right)={{s}_{1}G}_{1}\left(t;{t}_{1},\,{\sigma }_{1}\right)+{{s}_{2}G}_{2}\left(t;{t}_{2},\,{\sigma }_{2}\right)+{a}_{0}{e}^{-(t-{t}_{3})}$$
(5)
where C(t) is the concentration of fluorescein as a function of time t. The first two terms, Gaussian functions, represent the first and second pass of the bolus and the third term, an exponential decay, represents the washout phase51 The width of the first and second pass were σ1 = 10 s and σ2 = 25 s, respectively, and the decay rate of the washout phase, β, was 0.043 /minute. T1, t2 and τ are the time to peak for the first pass, second pass and washout phases, and were set at 0.171, 0.364 and 0.482 min, respectively51. S1, s2 and α were fixed at 0.833, 0.336 and 1.064 (dimensionless units). Peak concentration was normalised to unity at the inlet to the retinal artery and the time course in each connected vessel segment was time-shifted according to the velocity of blood in each vessel and scaled according to the ratio of flow in the parent and child vessels at bifurcation points.Image datasetsThis study was carried out in accordance with the Declaration of Helsinki71. Ethical approval of retrospective audit data was obtained through Moorfields Eye Hospital Research and Development and Audit number 1078. Informed consent was not obtained for retrospective data. Participants were not compensated as retrospective data was utilized. Only de-identified retrospective data was used for research, without the active involvement of patients. Sex was self-reported by study participants. Clinical ophthalmological retinal images were obtained from equipment at Moorfields Eye Hospital NHS Trust, London, UK: OCT-A images were obtained from a PLEX Elite 9000 (Carl Zeiss Meditec LLC, Dublin, CA, USA), ultra-wide true colour retinal photographs were obtained from Zeiss Clarus 500 Fundus machine (Carl Zeiss Meditec LLC, Dublin, CA, USA), fluorescein angiograms were obtained from Optos widefield camera (Optos, Inc. Marlborough, MA, USA). 19 manually segmented OCT-A images were obtained from healthy controls not ascertained for disease status (10 male, 9 female, mean age 39.89 (s.d. 11.25)). These manual segmentations were used in comparison of network structure with simulated networks. Datasets of 570 FA images, 590 colour retinal photographs, 43 OCT-A en-face images, and 130 simulated networks were used in training and testing the PI-GAN algorithm.Manual labelling of clinical dataManually labelled data was generated using a custom-built Python package enabling tracing of vasculature in 3D. The process involved placing user defining control points on the 2D image indicating where in a slab the vessel is located via maximum intensity projection. The z-height of the vessel was then fixed by identifying the height of the highest signal intensity voxel, which was manually constrained to exclude the choroid or RPE. The radius of each vessel was automatically calculated by setting a user-defined signal intensity threshold. Review of segmented structures was performed in 3D panel to assess and ensure labelling quality. In images with pathological blood vessels such as DR the abnormal vasculature or areas of neoangiogenesis were traced in the same manner. Vessel information (vessel coordinates, edge connectivity, number of edge points, edge point coordinates, radii, and vessel type) was exported and stored in Amira spatial graph format (ThermoFisher Scientific, Waltham, Massachussetts USA). Retinal regions were labelled. The macula was defined as a 5.5 mm diameter circular area centred on the fovea. The vessels surrounding the optic disc were labelled as a 3.6 mm diameter centred at the optic disc. Vessels outside these regions were defined as ‘peripheral’.Deep generative learningImage-to-image translation was performed using cycle-consistent generative adversarial networks30. This algorithm enables automated unsupervised training with unpaired samples, learning a bi-directional mapping function between two different domains with deep generative adversarial networks. It utilises cycle consistency, where the reconstructed image obtained by a cycle adaptation is expected to be identical to the original image for both generative networks. Cycle-consistent GANs are composed of two main deep neural network blocks which are trained simultaneously: an image generator (generator) and an adversarial network (discriminator). There is a loss (G loss) to make a synthesised image from domain A closer to a real image from domain B, and a loss (D loss) to distinguish the synthesised image from domain A from a real image from A. There are also losses facilitating the conversion in the opposite direction (G loss making synthesised image from domain B closer to domain A, and D loss to distinguish synthesised and real domain B images). Additionally, cycle loss is the difference between the input image and the double-synthesised image and identity loss is the difference between output and input images. A train/validation/test split of 75%/5%/20% was used. All PI-GANtraining and evaluation was performed using a single NVIDIA Titan RTX GPU.We iteratively trained a switchable PI-GAN algorithm with 500 epochs (Supplemental Fig. 5). All networks were trained using the optimizer ADAM solver35 with β1 = 0.5, β2 = 0.999. The learning rate for the first 100 epochs was 2 * 10−4, and then linearly decayed to 2 * 10−6. Images were pre-processed with crop size 256 pixels. The minibatch size was 1. The loss weights λ were set as 10. The model was trained on NVIDIA TITAN RTX in Pytorch v1.9.1.Loss functionsAdversarial lossesFor mapping function G : X → Y and its discriminator DY, the objective is expressed as:$${L}_{{GAN}}\left(G,\,{D}_{Y},\,X,\,Y\right)= {{\mathbb{E}}}_{y} \sim {p}_{{data}}\left(y\right)\left[log {D}_{Y}\left(y\right)\right]\\ +{{\mathbb{E}}}_{x} \sim {p}_{{data}}\left(x\right)\left[\log \left(1-{D}_{Y}\left(\right.G(x)\right)\right],$$
(6)
Where G attempts image generation G(x) to look similar to images from domain Y and DY tries to distinguish between translated samples G(x) and real samples y55.Cycle consistency lossCycle consistency loss is utilised to reduce the space of possible mapping functions to the target domain; the image translation cycle aims to bring an image back from its translation to the original. The images should satisfy both forward and backward cycle consistency. The mapping from one image, to translated domain, back to original is incentivised using the cycle consistency loss:$${{{{\mathcal{L}}}}}_{{cyc}}\left(G,\,F\right)=\, {{\mathbb{E}}}_{x \sim {p}_{{data}}\left(x\right)}\left[\left|\left|F\left(G\left(x\right)\right)-x\right|\right|1\right]\\+{{\mathbb{E}}}_{y \sim {p}_{{data}}(y)}\left[\left|\left|G\left(F\left(y\right)\right)-y\right|\right|1\right]$$
(7)
Identity lossAn additional loss is utilised to regularise the generator to be near an identity mapping when real samples from the target domain are used as the input of the generator:$${{{{\mathcal{L}}}}}_{{identity}}\left(G,\,F\right)={{\mathbb{E}}}_{y \sim {p}_{{data}}\left(y\right)}\left[\left|\left|G\left(y\right)-y\right|\right|1\right]+{{\mathbb{E}}}_{x \sim {p}_{{data}}\left(x\right)}\left[\left|\left|F\left(x\right)-x\right|\right|1\right]$$
(8)
Statistical evaluation of synthetic vessel networksVessel metrics of vessel branching angle, length, tortuosity, network volume and diameter were calculated. Two-sided Analysis of variance (ANOVA) was used to assess differences in these metrics by retina region (optic disc, macula, and periphery) and by status (healthy control and simulated network) (Supplementary Table 2) with eye (right OD/ left OS), participant sex, and OCT-A imaging scan pattern used as covariates.Expert evaluation of PI-GAN vesselsA clinical expert (consultant ophthalmologist at Moorfields Eye Hospital with fellowship-level subspecialty training in medical retinal disease and extensive clinical experience (14 years of experience)) reviewed the original segmentations from DRIVE and STARE datasets along with the corresponding images. Using evaluation in regions ‘disc’, ‘centre’ and ‘mid-periphery’, they evaluated all reviewed images as containing ‘possible’ small vessels. Their evaluation was that for 90% of DRIVE and 80.95% of STARE images, PI-GAN segmented additional vessels that were not seen in the original segmentation, that they considered ‘possibly’ present in the retinas. The remaining images were ranked ‘unclear’. They noted that image quality for DRIVE/STARE datasets was a limiting factor in conclusively determining the presence of smaller vessels.Evaluation metricsFrechet inception distanceGAN output was evaluated using the Fréchet Inception Distance (FID), which evaluates model quality by calculating the distance between feature vectors for real and generated images. FID compares the distribution of generated images with distribution of real images that were used to train the generator. Lower FID scores indicate more similarity between two groups. The FID score is calculated by first loading a pre-trained Inception v3 model. The output layer of the model is removed and the output is taken as the activations from the last pooling layer, a global spatial pooling layer.Three FID scores were calculated: real simulation images (synthetic space) versus manually segmented vasculature clinical images; real retinal photographs (retinal photograph space) versus PI-GAN generated retinal photographs; real OCT-A images (OCT-A space) versus PI-GAN generated OCT-A images; real FA (FA space) versus PI-GAN generated FA images.Dice scoreDice scores were additionally calculated. This is a commonly used performance statistic for evaluating the similarity of two samples. For a ground truth segmentation label L and associated prediction P, we measure the binary Dice score D:$$D\left(P,\,L\right)=\frac{2\left|L\cap P\right|}{\left|P\right|+\left|L\right|}$$
(9)
We carried out benchmarking of the PI-GAN algorithm against other models trained for manual segmentations from segmentation of retinal vessels using STARE, and DRIVE datasets public datasets, which are regularly used for benchmarking of algorithm results19,72. Dice score were evaluated from the output of PI-GAN trained to carry out the mapping between simulated data segmentations and retinal photographs and compared to GAN performance without synthetic data.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Physics-informed deep generative learning for quantitative assessment of the retina
