Learn Single-horizon Disease Evolution for Predictive Generation of Post-therapeutic Neovascular Age-related Macular Degeneration¹¹footnotemark: 1

Neovascular age-related macular degeneration (nAMD) is a main subtype of AMD. As the intravitreal vascular endothelial growth factor (VEGF) level elevates, choroidal neovasculars invade the avascular outer retina and severely damage photoreceptors, resulting in rapid vision loss [1]. At present, anti-VEGF injection is considered the preferred nAMD therapy [2]. Although ophthalmologists always give anti-VEGF injections after nAMD diagnosis, nAMD does not always respond satisfactorily to treatment. Besides, due to the lack of uniform guidelines, it is difficult for ophthalmologists to predict the short-term therapeutic response after anti-VEGF injection according to their subjective experiences [3]. This results in huge economic pressures and waste of resources [4]. Therefore, based on the known pre-therapeutic status of nAMD at time point $t_{1}$, predicting the post-therapeutic status of nAMD at time point $t_{2}=t_{1}+\Delta t$ can effectively forecast the efficacy of anti-VEGF injection for each patient. This can promote better clinical decision making.

Spectral-domain optical coherence tomography (SD-OCT) is a noninvasive, depth-resolved, high-resolution, and volumetric imaging technique. SD-OCT has become a pivotal diagnostic tool to visualize and quantitatively evaluate retinal morphological changes [5], including the diagnosis and tracing of nAMD [6, 7]. Each SD-OCT imaging can produce 3D volumetric images, also known as a cube. As shown in Fig. 1(a), we take Cirrus SD-OCT device as the example, each SD-OCT cube contains $1024\times 512\times 128$ voxels with a corresponding trim size of $2mm\times 6mm\times 6mm$ on the retina. In a cube, each slice along the vertical direction, with the size of $1024\times 512$, is known as a B-scan. A complete SD-OCT cube contains 128 continuous B-scans in space.

Diseases-associated predictions are more restrictive than general predictions. In the field of pattern recognition, multi-horizon predictions have been widely applied for natural language predictions, action predictions, video predictions, traffic predictions and etc. Given $\mathbb{X}_{{t_{1}}:{t_{N}}}=[{\bf X}_{t_{1}},{\bf X}_{t_{2}},\cdots,{\bf X}_{t_{N}}]\in{\mathbb{R}}$ as the historical N observations, each observation is obtained at the different time point and the time interval between any two adjacent time point is uniform. The actual future M observations are formally expressed as $\mathbb{X}_{{t_{N+1}}:{t_{N+M}}}=[{\bf X}_{t_{N+1}},{\bf X}_{t_{N+2}},\cdots,{\bf X}_{t_{N+M}}]\in{\mathbb{R}}$. We expect multi-horizon predictions to learn a mapping function $\mathcal{F}:\mathbb{X}_{{t_{1}}:{t_{N}}}\rightarrow\widetilde{\mathbb{X}}_{{t_{N+1}}:{t_{N+M}}}$ to obtain the prediction result $\widetilde{\mathbb{X}}_{{t_{N+1}}:{t_{N+M}}}$ as close as $\mathbb{X}_{{t_{N+1}}:{t_{N+M}}}$, as shown in Fig. 2(a). However, when applying multi-horizon predictions on medical images, several actual challenges raise:

• 1.

  Obtaining long-series observations from the same patient is intractable in practical clinical scenes, and it also is necessary to make effective predictions for a new patient with only one observation at the current time point.

• 2.

  For serial medical data, given any two different time points ${t_{i}},{t_{j}}\in({t_{2}},{t_{N}}]$, time intervals from adjacent time points may be different, namely ${t_{i}}-{t_{i-1}}\not\equiv{t_{j}}-{t_{j-1}}$.

• 3.

  Most medical observations generate 3D data and it is expensive for GPUs to learn from serial 3D data.

• 4.

  Therapeutic intervention dependent on medicine injection at a random time point is a key factor that cannot be ignored for diseases-associated predictions, and most medicine injection treatments only work for a short period of time.

• 5.

  It is well-known that the prediction accuracy decreases over time and risk tolerance on medical predictions is lower than other predictions, so clinicians always pay more confidence on the short horizon than the longer horizon.

These difficulties make it difficult for multi-horizon predictions to be really applied to medical images, thus learning a single-horizon prediction for medical images is more practical and realistic. Given one historical observation ${\bf X}_{t_{i}}\in{\mathbb{R}}$ at time point $t_{i}$, single-horizon prediction learns a mapping function $\mathcal{F}:{\bf X}_{t_{i}}\rightarrow\widetilde{{\bf X}}_{t_{i+1}}$ to obtain the prediction result at time point $t_{i+1}$, as shown in Fig. 2(b).

Nowadays, most existing disease prediction methods in the field of medical image processing fall into two classes, namely classification-based image-to-category (I2C) predictions, and regression-based image-to-parameter (I2P) predictions. Few works have focused on generation-based image-to-image (I2I) predictions, even though generative adversarial networks (GANs) have been widely used for modality transformation between medical images from different imaging devices. The aim of post-therapeutic prediction is to generate an image that explains how the anatomical appearance changes after treatment. We consider that the predictive post-therapeutic SD-OCT images would enable a better understanding of nAMD and clinical decision-making by presenting a visual post-therapeutic status of nAMD. According to the above, the target of our work is defined as:

Given a serial SD-OCT cubes with fixed time interval $\Delta t$, writing as $\mathbb{X}=[{\bf X}_{t_{1}},{\bf X}_{t_{2}},\cdots,{\bf X}_{t_{N+M}}]\in{\mathbb{R}}$, and any time point $t_{i}$ can be represented as ${t_{i}}={t_{1}}+(i-1)\Delta t$, $i\in[1,N+M]$. Anti-VEGF injection is given at each time point. For any SD-OCT cube ${\bf X}_{t_{i}}$, we hope to learn a mapping function $\mathcal{F}:{\bf X}_{t_{i}}\rightarrow\widetilde{\bf{X}}_{t_{i+1}}$ to predictively generate $\widetilde{\bf{X}}_{t_{i+1}}$ as close as ${\bf{X}}_{t_{i+1}}$.

However, learning the mapping function $\mathcal{F}$ directly based on a 3D SD-OCT cube is difficult and cost-consuming. Thus, for the SD-OCT cube ${\bf X}_{\it t_{i}}=[{\boldsymbol{x}}_{t_{i}}^{1},{\boldsymbol{x}}_{t_{i}}^{2},\cdots,{\boldsymbol{x}}_{t_{i}}^{128}]$, where $\boldsymbol{x}_{t_{i}}^{j}$ is j-th B-scan in ${\bf X}_{t_{i}}$, we learn the 2D mapping function $\mathcal{F}_{2D}:{\boldsymbol{x}}_{t_{i}}^{j}\rightarrow\widetilde{\boldsymbol{x}}_{{t}_{i+1}}^{j}$. $\mathcal{F}_{2D}$ is carried out 128 times until the SD-OCT cube is predicted completely.

In this paper, we present a Single-Horizon disease Evolution Network (SHENet) to predictively generate the post-therapeutic SD-OCT images based on pre-therapeutic SD-OCT images with nAMD. SHENet can help forecast the short-term response of anti-VEGF injection for individual nAMD patients. The main contributions in this paper are summarized as:

• 1.

  We explore the possibility of predictively generating post-therapeutic SD-OCT images based on pre-therapeutic SD-OCT images with nAMD, and further propose SHENet to solve this problem.

• 2.

  Graph evolution module (GEM) is proposed to imprison the process of disease evolution in the high-dimensional latent space by graph representation learning.

• 3.

  Evaluation reinforcement module (ERM) is proposed to reinforce the disease evolution process by combining an additional reconstruction generator and contrastive learning.

• 4.

  We design targeted experimental schemes according to actual clinical realities, and the results demonstrate that SHENet has great potential to generate post-therapeutic SD-OCT images with both high prediction performance and good image quality.

In clinical scenarios, several prediction requirements have been raised around AMD by ophthalmologists. In present, there are no uniform guidelines to be helpful for making predictions and ophthalmologists rely only on their own rich clinical experiences. Thus, ophthalmologists warrant the need for objective outcomes of subjective predictions. AMD typically develops from an early to an advanced form and advanced AMD is difficult to be cured effectively. When AMD is in its early stage, ophthalmologists predict the risk of progression to advanced AMD within the future short term in order to adapt therapies, recommendations, and follow-up frequency [8, 9, 10, 11, 12]. For advanced nAMD, ophthalmologists predict best-corrected visual acuity (BCVA) outcomes in patients receiving standard therapy [13, 14, 15]. For geographic atrophy (GA), as non-neovascular advanced AMD, there is a lack of effective treatments for retinal areas that have progressed to GA. But predicting the GA progression [16, 17, 18, 19, 20, 21] could allow for a better understanding of the pathogenesis and forewarn preventive treatment to normal retinal areas that are at high risk of developing GA in the future. Besides, several other predictions also reveal clinical needs. Bogunovic et al. [22] predicted low and high anti-VEGF injection requirements based on sets of SD-OCT images acquired during the initiation phase in nAMD. Liu et al. [23] and Lee et al. [24] exploringly generated individualized post-therapeutic SD-OCT images that could predict the short-term response of anti-VEGF injection for nAMD based on pre-therapeutic images using Pix2Pix. Forshaw et al. [25] predicted the visual gain from cataract surgery when the main cause of vision loss is nAMD. Pham et al. [26] generated future fundus images for early age-related macular degeneration based on generative adversarial networks.

Generative adversarial networks (GANs) [27] have become one of the widely leveraged techniques for generating images that look like the real thing. Traditional L1 or L2 supervision often results in blurred images [28], but GANs introduce an additional discriminator to play a min-max game with the generator that enforces the generator to output more realistic images. The discriminator is mainly used to estimate the divergence difference between generated fake images and real images. Different adversarial losses estimate different divergences between Wasserstein divergence distribution [29] and f-divergence family distribution [30]. If target domains own multiple distributions, conditional GAN (cGAN) [31] uses conditional labels to guide the generator to produce the images by fitting the specified distribution. For example, Yoo et al. [32] proposed a postoperative appearance prediction model for orbital decompression surgery for thyroid ophthalmopathy using a conditional GAN. GANs are also used for image-to-image (I2I) translation, in which pixel-level losses or feature-level losses are embedded to ensure the quality and stability of the generated image [26, 33, 34]. Besides, literatures [35, 36] explored the impact of discriminator architecture on GANs.

The random deviation of rotation angle and displacement among SD-OCT cubes at different time points caused by human operations cannot be learned, which further limits the I2I prediction. In other words, for two original SD-OCT cubes from the same patient captured at different time points, the B-scan with the same index in two cubes may not correspond to the same anatomical tissue. Thus, we perform voxel-wise serial alignment, including image flattening in vertical direction and image alignment in horizontal-axial directions, for all SD-OCT cubes before running SHENet. In this way, all SD-OCT cubes at different time points are aligned in the 3D space. The aligned SD-OCT cube can be seen in Fig. 1(b).

Image Flattening in Vertical Direction. We first obtain the locations of Bruch’s membrance (BM) of all B-scans by layer segmentation approach [37], and then all B-scans are flattened based on BM. Finally, we crop the B-scan restricted to the region from 0.75mm above BM to 0.25mm below BM. In this processing, negligible vitreum regions and sclera regions are removed as much as possible. Image flattening also reduces the size of the input SD-OCT images to alleviate the memory pressure of hardware.

Image Alignment on Horizontal-Axial Directions. We first generate the 2D vessel fundus images of each SD-OCT cube by restricting the projection region between inner segments/outer segments (IS/OS) and BM. The vessel fundus images are aligned using a scale-invariant feature transform (SIFT) flow method to obtain the transformation matrixes. Lastly, the transformation matrixes are applied to horizontal-axial directions of SD-OCT cubes to obtain the aligned SD-OCT cubes.

Pix2Pix [38] has been a great success to apply conditional GAN (cGAN) to supervised I2I translation. It regards input images as additional conditions to learn an I2I mapping, consequently producing specified output images. SHENet further extends Pix2Pix for our nAMD prediction task and its overview is illustrated in Fig. 3. In terms of model architecture, SHENet consists of:

(1) Prediction Generator $\mathcal{G}^{p}$ is the core of SHENet, including a feature encoder, a graph evolution module (GEM), and a feature decoder, which can predictively generate post-therapeutic SD-OCT images by inputting pre-therapeutic SD-OCT images with nAMD.

(2) Reconstruction Generator $\mathcal{G}^{r}$ is an auxiliary generator in the training process and will be removed in the model inference stage. $\mathcal{G}^{r}$ removes GEM from $\mathcal{G}^{p}$ and utilizes a reconstruction task to help $\mathcal{G}^{p}$ imprison the process of disease evolution in the high-dimensional latent space, and further to distill the function of feature encoder and decoder.

(3) Evolution Reinforcement Module (ERM) reinforces the process of disease evolution by working with $\mathcal{G}^{r}$ based on contrastive learning.

(4) Quality Discriminator $\mathcal{D}^{q}$ ensures that the predicted and reconstructed images look realistic.

(5) Pair Discriminator $\mathcal{D}^{p}$ determines whether predicted and reconstructed images are paired with input images in terms of pathological characterization that reflects therapy response.

The prediction generator $\mathcal{G}^{p}$ of SHENet consists of three modules: feature encoder ($\mathcal{G}_{E}$), graph evolution module (GEM), and feature decoder ($\mathcal{G}_{D}$). The pre-therapeutic SD-OCT images are mapped to high-dimensional latent space by $\mathcal{G}_{E}$, then GEM predicts the process of disease evolution after treatment, and finally, $\mathcal{G}_{D}$ recovers the predicted features to post-therapeutic SD-OCT images.

Feature Encoder ($\mathcal{G}_{E}$). Fig. 4(a) illustrates the architecture of feature encoder $\mathcal{G}_{E}$ that stacks several encoding blocks (EncBlock) and a mapping layer. Each EncBlock consists of two 3$\times$3 convolutional layers with ReLU activation, one channel attention block (CAB) [39], and one max-pooling layer for down-sampling. CAB (Fig. 4(b)) improves the quality of features by explicitly modeling the interdependencies between the channels of its convolutional features. Through CAB, informative features are selectively emphasized and less useful ones are suppressed.

Graph Evolution Module (GEM). Given the encoding features ${\boldsymbol{f}}_{t_{1}}^{j}\in{\mathbb{R}}^{{h^{\prime}}\times{w^{\prime}}\times{d^{\prime}}}$ by $\mathcal{G}_{E}$, where $P={h^{\prime}}\times{w^{\prime}}$ is feature numbers and ${d^{\prime}}$ is feature dimension. We consider the evolution of each feature should be related to other features. Intuitively, a series of convolution layers or a linear regression is easier to implement to predict the change in the feature level. However, a series of convolution layers or a linear regression defaults that all other features contribute equally to the targeted feature. In fact, the contributions of all features are unequal due to their different semantic information. Graph neural networks can model the contributions of all features automatically, which is more reasonable than a series of convolution layers or a linear regression. Let each feature be the vertex of the graph, and we can represent ${\boldsymbol{f}}_{t_{1}}^{j}$ as a fully-connected undirected graph $\mathcal{R}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ is vertex set and $\mathcal{E}$ is edge set. Further given the adjacency matrix ${\boldsymbol{A}}$, the diagonal degree matrix ${\boldsymbol{D}}$ and the identity matrix ${\boldsymbol{I}}$, the relation of features can be extracted by the following formula:

where ${\boldsymbol{W}}^{l}\in{\mathbb{R}}^{{d_{in}}\times{d_{out}}}$ is the learnable weight, and $\sigma$ is the Mish function [40]. ${\boldsymbol{H}}^{l}\in{\mathbb{R}}^{{P}\times{d_{in}}}$, ${\boldsymbol{H}}^{l+1}\in{\mathbb{R}}^{{P}\times{d_{out}}}$ are the input features and the updated features at l-th layer. Intuitively, for p-th vertex $\mathcal{V}_{p}$, its all neighbors contribute unequally to the evolution of $\mathcal{V}_{p}$. To model this, we introduce the multi-head GAT [41] to explicitly consider the importance of the neighbors. Take independent GAT as an example, we calculate the attentive score ${\gamma}_{pq}$ of the vertex pair $(\boldsymbol{h}_{p},\boldsymbol{h}_{q})$:

where $\mathcal{N}_{p}$ is the neighbors of p-th vertex in the graph, and $[\cdot]$ is a concatenation operation. $\boldsymbol{W}\in{\mathbb{R}}^{{d}_{out}\times{d}_{in}}$ indicates a linear mapping and $\boldsymbol{\alpha}\in{\mathbb{R}}^{{2d_{out}}}$ indicates a single-layer fully-connected layer. LReLU is the nonlinear activation. Then, we compute the average of multi-head GAT for the output of vertex $\mathcal{V}_{p}$:

where G=5 is the number of heads. Similarly, we repeat the GAT computation on each vertex to obtain the complete output. Considering the powerful ability of graph neural networks for information inference, GEM is only composed of 3 multi-head GATs with the channel numbers 1024, 1024, 1024. The overview of GEM is shown in Fig. 4(c).

Feature Decoder ($\mathcal{G}_{D}$). Fig. 4(d) illustrates the architecture of feature decoder $\mathcal{G}_{D}$ that stacks several decoding blocks (DecBlock) and a mapping layer. Each DecBlock consists of one de-convolution layer with ReLU activation, and two convolutional layers with ReLU activation.

To reinforce the process of disease evolution, we introduce reconstruction generator $\bf{\mathcal{G}}^{r}$ that removes GEM from $\bf{\mathcal{G}}^{p}$ to reconstruct $\overline{\boldsymbol{x}}_{t_{2}}^{j}$ by inputting ${\bf X}_{t_{2}}^{[j;\Delta s]}$:

Prediction generator $\bf{\mathcal{G}}^{p}$ and reconstruction generator $\bf{\mathcal{G}}^{r}$ share feature encoder ${\mathcal{G}}_{E}$ and feature decoder ${\mathcal{G}}_{D}$. We use the project head ${\it h}(\cdot)$ [42] to map ${\boldsymbol{f}}_{t_{1}}^{j}$, $\widetilde{\boldsymbol{f}}_{t_{2}}^{j}$, ${\boldsymbol{f}}_{t_{2}}^{j}$ to ${\boldsymbol{z}}_{t_{1}}^{j}$, $\widetilde{\boldsymbol{z}}_{t_{2}}^{j}$, ${\boldsymbol{z}}_{t_{2}}^{j}$ for evolution reinforcement learning:

where $\boldsymbol{W}^{(1)}$, $\boldsymbol{W}^{(2)}$ are weight matrixes, $GAP(\cdot)$ indicates global average pooling. We hope predicted $\widetilde{\boldsymbol{z}}_{t_{2}}^{j}$ and real ${\boldsymbol{z}}_{t_{2}}^{j}$ should be as similar as possible, while predicted $\widetilde{\boldsymbol{z}}_{t_{2}}^{j}$ and real ${\boldsymbol{z}}_{t_{1}}^{j}$ should be different. Thus, we build the evolution reinforcement learning based on contrastive loss:

where $\tau$=1 is the temperature factor, $Sim(\cdot)$ is the cosine similarity metric. ERM also further distills the function of the feature encoder and feature decoder.

For the predictively generated $\widetilde{\boldsymbol{x}}_{t_{2}}^{j}$ and the reconstructed $\overline{\boldsymbol{x}}_{t_{2}}^{j}$, we firstly use a quality discriminator $\mathcal{D}^{q}$ to ensure their image quality is realistic:

Furthermore, we consider that pathological manifestation of $\widetilde{\boldsymbol{x}}_{t_{2}}^{j}$ and $\overline{\boldsymbol{x}}_{t_{2}}^{j}$ should show the response of $\boldsymbol{x}_{t_{1}}^{j}$ after anti-VEGF injection. In other words, $\widetilde{\boldsymbol{x}}_{t_{2}}^{j}$ and $\overline{\boldsymbol{x}}_{t_{2}}^{j}$ should be paired with $\boldsymbol{x}_{t_{1}}^{j}$. Thus we also use a pair discriminator $\mathcal{D}^{p}$ to make the decision:

In SHENet, both $\mathcal{D}^{q}$ and $\mathcal{D}^{p}$ follow the popular PatchGAN [38].

In the training process, we use $\mathcal{L}_{L1}$ measures the L1 distance between output images and ground truth:

A pair discriminator $\mathcal{D}^{p}$ makes the decision whether ($\widetilde{\boldsymbol{x}}_{t_{2}}^{j}$, $\boldsymbol{x}_{t_{1}}^{j}$) and ($\overline{\boldsymbol{x}}_{t_{2}}^{j}$, $\boldsymbol{x}_{t_{1}}^{j}$) are paired:

A quality discriminator $\mathcal{D}^{q}$ ensures $\widetilde{\boldsymbol{x}}_{t_{2}}^{j}$ and $\overline{\boldsymbol{x}}_{t_{2}}^{j}$ are realistic:

We combine all losses together and the complete optimization can be represented as:

According to clinical research, the therapeutic effect of anti-VEGF injection for nAMD lasts about one month, so nAMD patients treated should be followed up one month later to evaluate the therapeutic response of nAMD. Then, ophthalmologists can make the decision on whether further treatment is needed. However, to the best of our knowledge, there are no large-scale, public, and annotated SD-OCT datasets that can be acquired for our model validation due to the differences in imaging protocols, privacy problems, and lack of medical integration. This is also a common problem in the field of medical image processing.

The experimental data includes 46208 paired SD-OCT images obtained from 383 SD-OCT cubes of 22 nAMD patients. Only one eye per nAMD patient is included in the dataset, namely 22 eyes are included. Each patient contains about 17 serial SD-OCT cubes captured at different time points. During the SD-OCT imaging performed once a month, ophthalmologists gave anti-VEGF injections for nAMD treatment. The time interval between any adjacent time points is about one month, that is to say, the predictive single horizon is fixed to one month. These SD-OCT cubes were captured by Cirrus SD-OCT device Carl Zeiss Meditec, Inc., Dublin, CA. The size of an SD-OCT cube is $1024\times 512\times 128$. Detailed attribute descriptions are presented in Table 1.

For the predictively generated post-therapeutic SD-OCT images, qualitative evaluation using human subjective judgment is the most straightforward and effective way to verify the prediction performance of SHENet. Besides, we also use three metrics, peak signal-to-noise ratio (PSNR), structural similarity (SSIM) and learned perceptual image patch similarity (LPIPS), to evaluate the prediction performance quantitatively. These three metrics have been previously used for the evaluation of image generation [43]. PSNR measures the image quality and higher PSNR means less distortion. SSIM measures the structure similarity and higher SSIM means higher structure similarity. LPIPS measures the similarity of deep visual features and lower LPIPS means higher feature similarity. To make LPIPS keep the consistency with PSNR and SSIM, namely higher value indicates better performance, here we use 1-LPIPLS to replace LPIPS.

To highlight the advantages of SHENet in terms of I2I prediction performance, we choose four GANs-based methods for competitive comparison, namely Pix2Pix [38], Pix2PixHD [44], Fundus2Angio [45], Att2Angio [46]. Pix2Pix and Pix2PixHD have become popular baselines and are widely used for image generation. Fundus2Angio and Att2Angio are the latest GANs-based methods for modality transformation in the field of medical image processing.

In order to conform to real clinical application scenes, we design our experiments from three sub-evaluations:

P-0 Evaluation. For the new nAMD patients that only have one SD-OCT imaging at time point $t_{1}$, we want to predict post-therapeutic SD-OCT images at time point $t_{2}$. Thus, the SD-OCT images from all other patients are used for training SHENet. We use five-fold cross-validation until all patients are tested, as shown in Fig. 5(a).

P-1 Evaluation. For the nAMD patients that have two SD-OCT imaging at time points $t_{1}$ and $t_{2}$, we would like to predict post-therapeutic SD-OCT images at time point $t_{3}$. Thus, we transfer the model parameters from P-0 evaluation and fine-tune SHENet only using $({\bf X}_{t_{1}},{\bf X}_{t_{2}})$. In the model inference stage, we take ${\bf X}_{t_{2}}$ as model input to predictively generate ${\bf X}_{t_{3}}$. We use five-fold cross-validation until all patients are tested, as shown in Fig. 5(b).

P-M Evaluation. For the nAMD patients that have owned multiple regular SD-OCT imaging, we would like to predict consequent post-therapeutic SD-OCT images. Thus, we remain the last two SD-OCT cubes of all patients for model validation and the remaining SD-OCT cubes are used for training model, as shown in Fig. 5(c).

The experiment is constructed in a hardware condition with Intel Xeon CPU, one GeForce RTX 3090 GPU and 128 GB RAM, and a software condition with Python3.5 and Pytorch.

For the input multiple B-scans, we choose $\Delta s=3$ and zero-padding operation is used if $j-\Delta s<0$ or $j+\Delta s>128$, thus the model input is a three-channel SD-OCT image. The number of encoding blocks and decoding blocks is 5. The output dimensions of 5 encoding blocks are {128, 256, 512, 1024, 2048} respectively. The output feature size of the feature encoder is $16\times 16\times 2048$, which means the vertex number of a fully-connected graph in GEM is $16\times 16$. In the loss function, $\lambda=100$ and $\mu=10$ balance the weights among GAN losses, L1 losses and contrastive loss. Flip and rotation operations are used here for data augmentation.

Adam optimizer with initial learning rate of 0.0001 and weight decay of 0.1 is chosen for model optimization. The batch-size is set to 2. SHENet was trained for 100 epochs and the training was properly completed.

We qualitatively compare our SHENet with competing methods by visualizing the examples from two different nAMD patients based on P-0, P-1 and P-M evaluations, as shown in Fig. 6. In terms of image quality, benefiting from adversarial training, all methods can produce clear and unblurred SD-OCT images. Besides, SHENet can stably maintain the structural integrity of the retina to achieve the better visual effects and actually predict the status of nAMD one month later after anti-VEGF injections.

We quantitatively evaluate the prediction performance based on P-0, P-1 and P-M evaluations using PSNR, SSIM and 1-LPIPS metrics, as shown in Table 2. For P-0 and P-1 evaluations, we calculate the average of five-fold cross-validation as the final results. Overall, SHENet obtains consistently superior results than competing methods on three experimental designs. Among three experimental designs, SHENet achieves the best prediction performance on P-M evaluation (PSNR: 24.198, SSIM: 0.349, 1-LPIPS: 0.642), followed by P-1 evaluation (PSNR: 23.875, SSIM: 0.337, 1-LPIPS: 0.626), and P-0 evaluation (PSNR: 23.659, SSIM: 0.326, 1-LPIPS: 0.609) has the worst results, which shows the consistency with qualitative evaluation in Fig. 6. That is because the difference among patients is significantly greater than that among serial SD-OCT observations from the same patient. We consider that the gap will narrow by further collecting more samples from more patients.

Evaluation of Hyper-parameter Values. We analyze the influence of hyper-parameters $\lambda$ and $\mu$ in loss function with the experiment. We first fix $\mu$ to be 1 and vary $\lambda$ from 10 to 190 with interval 10. As shown in Fig. 7(a), 1-LPIPS metrics of three evaluations are consistently increasing when rising the of $\lambda$ from 10 to 100 and consistently decreasing when further rising $\lambda$. We then fix the value of $\lambda$ to 100 and vary $\mu$ from 1 to 19 with interval 1. As shown in Fig. 7(b), the performance of SHENet achieves the highest values when $\mu=10$.

Evaluation of Model Inputs. In terms of model input, we investigate the difference between single B-scan and multiple B-scans, and the quantitative comparisons are recorded in Table 3. We find that the prediction results using multiple B-scans as model input demonstrate significant improvements to those using a single B-scan as model input on three experimental designs.

Evaluation of Model Architecture. In the training process, we build our single-horizon disease evolution in the high-dimensional latent space based on the cooperation of GEM, $\mathcal{G}^{r}$ and ERM. We investigate the impact of three components by stacking them one by one. As shown in Table 4, SHENet achieves better results than others when considering these three components jointly. This evidences that SHENet effectively imprisons the process of disease evolution in GEM and further reinforces the process by $\mathcal{G}^{r}$+ERM. Therefore, adopting all three can improve prediction performance.

Evaluation of Discriminator. To verify the relevance of two discriminators ($\mathcal{D}^{q}$, $\mathcal{D}^{p}$), we separately analyze the prediction results of SHENet when one of them is reserved. When only using one of the two, quantitative metrics in Table 5 show performance degradation. Besides, we also observe that single $\mathcal{D}^{p}$ performs better than single $\mathcal{D}^{q}$, as $\mathcal{D}^{p}$ can also control the image quality to a certain extent. Thus, it shows the necessity to use two independent discriminators to respectively control the image quality and the pathological characterization.