Photonic Reservoir Computing for Real-Time Respiratory Motion Prediction

: The application of reservoir computing (RC) is for the first time studied in a class of forecasting tasks in which signals are under random physical perturbations, meaning that the data-baring waveform distortions are versatile, and the process is not repeatable. Tumor movement caused by respiratory motion is such a problem and real-time prediction of tumor motion is required by the clinical radiotherapy. In this work, a true-time delay (TTD) respiration monitor based on photonic RC with adjustable nodes connection is developed specifically for this task. A breathing data set from a total of 76 patients with breathing speeds ranging from 3 to 20 breath per minute (BPM) are studied. A double-sliding window technology is demonstrated to enable the real-time establishment of an individually trained model for each patient and the real-time processing of live-streamed tumor position data. Motion prediction of look-ahead times of 66.6 ms, 166.6 ms and 333 ms are investigated. With a 333 ms look-ahead time, the real-time RC model achieves an average normalized mean square error (NMSE) of 0.0246, an average mean absolute error (MAE) of 0.338 mm, an average therapeutic beam exposure efficiency of 94.14% for an absolute error (AE) < 1mm and 99.89% for AE < 3mm. This study demonstrates that real-time RC is an efficient computing framework for high precision respiratory motion prediction.

Reservoir computer (RC), often referred to as an echo-state network 1 , is one class of recurrent neural networks (RNN) that is rather efficient at processing temporal or sequential signals. The neurons (nodes) in the hidden layer of a RC have random but fixed connections, and only the readout layer weights are trained using computationally inexpensive linear regression methods 2,3 . Compared to machine learning frameworks of interactive learning, a RC is an appealing artificial neural network (ANN) approach for its simplicity in network architecture and algorithm implementation, and thus fast computation capability.
To date, RC has been applied to both real life and benchmark computation tasks such as pattern classification 4 ,5 and generation 6 , time series forecasting 7 , voice recognition 8 , equalization for wi r el ess channel 9,10,11 , satellite communication 12 and so on. Recently, a photonic RC has been put to the test in a classification task, more specifically optical signal recovery 13 in fiber communication that requires ultrafast signal processing speed to classify largely distorted optical signals. In this task, the signal waveform distortion is caused by multiple sources of nonlinear perturbations when the optical signals travel along the fiber. A common theme here is that the perturbations, though intense and nonlinear, are applied constantly to the traveling signals so the waveform distortion is repeatable when the signals go through the same fiber channel. There exists another class of problems, little studied so far, in which the physical perturbations to the signals are stochastic, resulting in non-repeatable waveforms. Living biological systems often behave in this manner and can display greatly varying amplitudes at random. Extreme signal distortion diversity and randomness complicate the ability of RC to recognize and classify signal patterns.
In the medical field, respiration induced inner organ/tumor motion provides an example of a random temporal signals where the tumor follows the patient's breathing pattern, but in a complex, quasi-sinusoidal motion. The breathing patterns are intermittent as they are subject to both the natural variation of the patient's respiration and the influence of the patient's physical conditions at any given moment, such as slow upper body movement, compromised lung function and pain or nervousness. The tumor's respiratory motion can also exhibit extreme irregularities during involuntary events such as sneezing and coughing, causing large displacement variations in amplitude and frequency from the baseline quasi-sinusoidal respiration curve, and sometimes even inducing spikes in the temporal breathing pattern. In a clinical deployment, the respiratory motion data is collected along with the activated radiation beam. Tumor movement requires that the therapeutic beam be adjusted in vertical and lateral position so as to always remain on-target. This tumor tracking method is commonly referred to as adaptive radiotherapy. Several machine learning methods have been explored for respiratory motion prediction 14 , 15 , 16 , 17 and have demonstrated effectiveness. However, none of these learning-based approaches can delivery satisfactory prediction efficiency with the required measurement latency. The forward neural network (FNN) approach does not account for residual memory effects from previous inputs during the training epoch. On the other hand, RNNs can process temporal signals with feedback from earlier inputs 18 , 19 , producing excellent prediction precision. The drawback of classic RNN is that they are notoriously difficult to train. For example, it can take over 20 hours to develop a training model using a long-short term memory (LSTM) algorithm for motion prediction 20 .
The respiratory motion signals are represented by nonstationary data and forecasting the tumor motion trajectory in a real-time frame is essential for clinical practice. In this work, we present a photonic RC approach that provides a valid solution for real-time respiratory motion prediction. The reservoir layer of the RC is implemented in a true-time delay (TTD) photonic circuit constructed from off-the-shelf commercial components. A RC algorithm and strategy that lead to the least computation expense is of paramount importance for real-time motion prediction. In this work, we studied a double-sliding window data processing strategy in which the input data set window and the training data set window slide together with sufficient intervening time between the two operations to process the livestreamed tumor position data without discernable latency. The training data set has only 600 data points for fast algorithms execution while sliding the training data set ensures that the weight matrix is trained only on the latest respiration signals while discarding distant ones. An enhanced precision in the predicted position of the tumor is thus obtained in a real-time frame as compared with predictions obtained with a RC network trained on a fixed data set. Unsynchronized TTD photonic RC is implemented in order to increase the dynamics of the reservoir and reduce the reservoir state vector dimension in favor of fast computing. In this work we show that tumor motion can be precisely predicted using an algorithm based on a photonic true-time delay reservoir (TDR) computer. Furthermore, it is shown that such a computer can yield a generalized processing approach for a class of nonstationary temporal signals characterized by stochastic perturbations and a r ich diversity waveform distortion.

Delay-line RC network
Several RC topologies have been reported based on how the neurons in the hidden layer are connected 21,25 . Here, we implement a hidden layer delay line topology with adjustable connections between nodes in a fiber-optic TTD implementation. A single nonlinear physical node is utilized to map and project the input layer signals to the reservoir layer. A large number of virtual nodes are stored temporally along a delay line. A schematic representation of the reservoir computer is shown in Fig. 1. A position sensor records and sends digitized sequential motion signals every = 33.3ms. The digitized input signal ( ) undergoes a sample-and-hold operation in a duration marked by ' to produce a continuous function ( ), where n is the motion data index. The input signal ( ) is then projected to the reservoir layer by multiplying with a mask function ( ) with a periodicity identical to ′. The piecewise mask function levels are drawn from a random distribution such that the duration of each piece, , and the number of virtual nodes, ', satisfies the relation = ′ / ′ . The photonic RC processes continuous analog signals that have the form of ( ) = ( ), ′ ≤ < ( + 1) ′ (1) where represents the ℎ node in the delay line. The mapping of low dimensional input data to a reservoir state in a high dimensional dynamic space greatly enhances the separability of classes of data 22 . The reservoir state values ( ) are impacted by both the current inputs and the past inputs, and can be described by an evolution equation: where is the reservoir internal connectivity; and α, β and ϕ are the hyperparameters of the RC that needs to be optimized for one particular computation task. The delay line in the reservoir forms a closed loop, similar to a feedback loop or an oscillator cavity 23 . The round-trip delay, i.e., the oscillator characteristic time, denoted by , is determined by the combined optical and electrical signal traveling time in the loop. The reservoir state can be discretized by sampling a single point on the delay line, with a period , to produce ' pieces of nonlinear response, ( ), to each input data ( ). The dynamics of the reservoir for a typical synchronized RC operation 5 , i.e., ′ = , can be described by discretized reservoir states ( ): An optical Mach-Zehnder modulator (MZM) is used as the nonlinear node and provides the sinusoidal transfer function. Mathematically, other nonlinear functions, such as hyperbolic tangent function can also be used 24 for the RC. The projected signals in the reservoir travel unidirectionally in the delay-line loop and the optical power attenuated until the signal strength is reduced below the noise level, emulating short-term memory of a biological neural system, which will be referred to in this work as the decaying memory. The target output, i.e., the ground truth, is marked as ( ) and �( ) is the predictive value. The output weight matrix is calculated via ridge regression 25 during training: where the subscript T denotes the matrix transpose operation, is the regularization factor typically very small, is the target output matrix, is the reservoir state matrix and is the identity matrix. The most computationally expensive step is matrix inversion operation. Once is determined in the training process, the RC response �( ) is the linear sum of the weighted reservoir output neurons: Fig. 1: A schematic representation of photonic a TTD topology. A total of virtual nodes are connected in a closed-loop ring architecture. The input signals are projected to the reservoir via a mask function while a gain is set for the coupling strength. The nonlinear node, Sin, provides a sinusoidal function with an internal loop signal amplification factor (gain) of .

RC algorithms for real-time tumor position prediction
Current machine learning methods often require acquisition of a minute to several minutes of motion signals to properly train a neural network for the required computational precision 20,26,27 . Clinic practice requires a maximum desirable training time of less than 30 seconds. This is because the radiation beam needs to be shuttered during neural network training and beam re-alignment needs to be verified. This adds overhead time to the total radiation treatment efficiency. In this work, we trimmed the training data set length so as to decrease the training data length from several minutes to 20 seconds with improved predictive location precision in a real-time clinical setting. The training data length ( ) is defined as the number of data points used for RC network training. A smaller not only enhances the radiation treatment efficiency but also reduces the reservoir state dimension favoring faster RC algorithm execution.
A novel dual-sliding window strategy is adopted in the execution of the algorithm, as exemplified in a small portion of a respiratory motion curve shown in Fig. 2. The training data set , having fixed length, slides along the breathing curve. This generates a new training data set with every newly acquired position sensor data, thus establishing an updated matrix with time. This strategy enables uninterrupted, real time, respiratory motion prediction with minimal therapeutic beam shuttering, resulting optimal beam utilization time and optimal clinical beam dose efficiency.
We evaluated between 300 and 1000 data points and observed that the = 600, corresponding to = 20 seconds, is sufficient to produce satisfactory motion prediction results for most breathing patterns. The second sliding window refers to a fixed input data length, used in forecasting tumor position, is defined as the input data length . A larger means more position data points from the past are utilized in forecasting of the next several tumor positions. Typically, a greater value produces better prediction results; however, for > 15, the increase in prediction precision is diminished so is kept at 15 in this work. The look-ahead time, is defined as a future time point at which the respiration caused tumor movement is calculated by the RC. The look-ahead time ultimately is determined by the electro-mechanical response time of the leaf adjustment mechanism used to define the beam aperture (see Supplementary for details). The number of data points within the look-ahead time duration is ℎ . In this work, we explored three look-ahead times: = 66.6 ms, 166.6 ms and 333 ms corresponding to ℎ = 2, 5 and 10, respectively. Fig. 2: A schematic representation of the algorithm implemented in our photonic RC network for real-time tumor motion prediction. A dual-sliding window comprising (~20 sec) and (=15 data points) is used in this figure. The maximum look-ahead time studied in this work is = 333 ms which corresponds to ℎ = 10. The 15 data points in green dots represent the data of the input window that are used to predict the data point circled in red 333ms into the future. The red dot represents the data to be forecasted in a future time of 333 ms.

Results
A total of 76 breath patterns were collected and analyzed for motion prediction using the photonic TDR. Signal noise in data acquisition and digitization is inevitable due to the electronic components used in the devices including the position sensor and other electronic instruments. A digital Gaussian filter is applied on the motion data to minimize these noise effect (see more details in supplementary material).
Respiratory motion patterns of humans are inherently complex and irregular. They vary with time, differ drastically among individuals, and often include abrupt changes due to events such as coughing and sneezing. It is also common that a breathing pattern of a patient combines multiple irregularities, resulting in exceedingly versatile temporal motion curves. Mathematically, temporal signals with a large degree of waveform distortion from a baseline quasi-sinusoidal function, would require much higher class separability in the high-dimensional reservoir space in order to achieve the same level of prediction precision. For adults at rest, the respiration rate is normally at 12-16 breaths per minute (BPM). In a full breath cycle, the patient goes through inhale and exhale phases, leading to a motion curve of peak-valleypeak. At cancer radiation clinics, a wide range of BPM values, ranging from 3-20 28,29 is often observed and treated by radiotherapy. The 76 breathing patterns (BP) are arranged in a sequence from low to high BPM in Fig. 3. The prediction errors are evaluated by the normalized mean square error (NMSE) and the mean absolute error (MAE), defined as where � stands for the predictive value and � is the averaged position value. As the position sensor collects data at a fixed rate of 30 Hz, a higher BPM implies a lower density of data points per respiration cycle. The NMSE and MAE variations with BPM are plotted in Fig. 3 for all 76 breathing patterns in a real-time prediction frame. Both breathing irregularities and BMP affect prediction precision. A trend line to the 5 th order polynomial is also plotted in Fig. 3. A gradual increase of prediction error with BPM for both NMSE and MAE is observed as a higher BPM is equivalent to a smaller sampling rate of the breathing waveform. Therefore, a simple solution to improve the overall prediction precision is to increase the position sensor sampling rate. For BMP > 18, the trend line bends downwards slightly. This can be explained as the randomness in the waveform distortion. Only a few patient breathing curves fall within the range outside a BMP > 18 so there also exists statistical uncertainty. We have selected 15 representative motion curves with different irregularities and plotted the predictive values in BPM groups of 1-10, 10-15 and 15-18, as described in the caption to Fig. 4. All three look-ahead times are evaluated for the NMSE and MAE. In one respiration cycle, the largest error occurs near a peak or valley. Increased prediction errors are also observed near breathing irregularities that cause localized large motion displacement such as spikes. Absolute errors (AE), i.e., = | − �| are also plotted in Fig.  4 to monitor the real-time AE. Clinically, an AE < 1 mm for the beam exposure margin is considered outstanding, AE < 3 mm is adequate and AE < 5mm is acceptable 30 . Three horizontal lines are plotted to mark the 1 mm, 3 mm and 5 mm margins in the AE plot.

Discussion
RC for improved radiation delivery efficiency. In a clinical setting, the real-time photonic TDR respiratory motion prediction strategy and instrumentation described here is anticipated to be adopted in a dynamic tracking system used in combination with the respiratory gating for radiotherapy 31,32,33 . Radiation is synchronized with the gating window and the beam is only de-shuttered when the position prediction error is lower than a pre-set error margin. The ratio of the gating window over one respiration cycle marks the radiation treatment efficiency. NMSE, MAE, or RMSE (root mean square error) 34 are commonly used figures of merits of statistics. In this case, these figures of merits are used to assess machine learning outcomes while they are not explicitly used for evaluating the radiation treatment efficiency. In this work, the therapeutic beam efficiency (TBE) via the time percentage of the predictive value with AE < 1 mm, 3 mm, 5 mm margins within one radiotherapy session is calculated and plotted in Fig. 5. For a majority of the 76 motion curves used in this study, a therapeutic beam exposure efficiency of 80% can be obtained for AE < 1 mm for all three look-ahead times while > 98% efficiency for AE < 3 mm margin. In comparison, the current clinically implemented respiration gating method is typical to set the gating phase of 20% to 70% 33 . Assuming equal time in each phase, the radiation beam is on for phases from 0% to 20% and 70% to 90%, giving ~ 50% as the radiotherapy beam exposure time efficiency. Using the real-time RC algorithm implemented on a photonic TDR, a significant increase in beam delivery efficiency is observed for both small and large AE margins in all three look-ahead time conditions.

Real-time RC versus fixed training model prediction.
To date, learning-based respiration motion prediction algorithms often use a pre-collected breathing pattern library to train the neural network to lift the constraint of in-line training due to large computational resource requirements. This training strategy is known to be subject to relatively larger errors as the breathing pattern characteristics tend to vary drastically among individuals or even for the same patient. This means that the respiratory motion training and testing data collected at different times is largely inadequate. Another approach is to use a portion of one breathing motion data set for training while the remaining data set is used for validation and prediction. In this case the training data yields a fixed that is used for all motion signal predication calculation. This approach is termed "fixed-model RC" in this work. For a real-time RC, varies as the training data window k tr and slides along the respiratory motion data. To compare the fixed-model RC with the real-time RC, the NMSE and MAE are calculated and plotted in Fig. 6. For the real-time RC, the sliding training signal window is selected to be = 20 sec ( = 600). For the fixed-model RC, two scenarios of = 20 sec and 60 sec ( = 1800) are used. For the latter case, 17 breathing motion data sets having a total signal length < 1 minute are excluded from the fixed-model RC calculation. TBE measures the radiation toxicity caused by the tumor respiratory movement, so an outstanding TBE is important for clinic deployment of the real-time RC. TBE for absolute error (AE) < 1mm, AE < 3mm and AE < 5mm are plotted in Fig. 7 for comparison. TBE in real-time has outperformed fixed model RC for all AE scenarios in both testing cases. For a tighter AE threshold, i.e., AE < 1mm, the improved prediction accuracy has drastically increased the TBE. For example, in Fig. 7(a), only 4 out of the 76 patents have TEB < 80% using real-time RC while 44 patients have TBE < 80% for the fixed RC model. Real-time RC has demonstrated better prediction precision in NMSE and MAE than the fixed-model RC approach for both scenarios, even though for the latter case, a longer training data set = 60 sec is used for the fixed-model RC. As the distortion mechanism of the waveform formed by the peak-valley-peak pattern of the respiratory motion is highly versatile, the breathing pattern characteristics vary continuously with time but can transition rapidly to a different pattern just in a few respiratory cycles. When distortions or irregularities appear frequently in the breathing waveform, earlier motion signals become less correlated to the current breathing pattern. An enhanced class separability, reflected by the computed Wout matrix can give more accurate predictive values by discarding those earlier data in the training. By updating the training model with the latest acquired position data, the real-time RC has the capability to track the pattern changes. The real-time RC also demonstrates considerably improved AE results than that of fixed-model RC, indicated by the TBE calculation.
Dynamics of delay-line reservoir. The separation property of a RC can often be improved by increasing the reservoir dynamics. It is commonly perceived that a larger number of virtual nodes in a RC will result in richer dynamics so improved class separability. Recently, a lemma from information theory that supposes the existence of an optimal reservoir node number for a given computation task, is analyzed mathematically in 27 . In the photonic TDR implementation, the virtual node number is adjusted by changing the fiber delay length, i.e., . Without any extra fiber, the photonic TDR has an intrinsic delay time of 28 ns, corresponding to = 14 in our setting. In the experiment, we studied the inclusion of an extra fiber up to 1km in length for extended true-time delay. A 1km fiber corresponds to a signal delay of ~ 4. 9 and = 619. The tested NMSE varies slightly in the range of = 14 to 619. As empirically the test is stable at small number of nodes, we selected = 14 with = 28ns for fast RC computing.
In a physical RC, the slowest component in the reservoir oscillator determines the response time (denoted as TD) of a TDR. The feedback photodetector (PD) has a marked bandwidth of 1 GHz, so the response time of the TDD is comparable to (2ns). The reservoir states of adjacent virtual nodes in the time domain would then interact in the form of electromagnetic wave resulting in complex transient response of the reservoir. The quantitative relation of the response time with was discussed in Appeltant et. al. 22 and Larger et. al. 35 . Unsynchronized virtual nodes can also lead to enriched reservoir dynamics 5 . In this work, the node separation is fixed at 2ns, choosing ' = 12 resulting in a node mismatch of = -′ = 2 for unsynchronization. Other node mismatch conditions and how affects the reservoir dynamics are described in detailed in the supplementary material.

Computation time consideration. For the real-time RC,
is computed and refreshed with any new motion data every = 33.3ms which sets the time budget for each computing operation. In general, the input signals are sent to the RC at an input rate set by an analog-to-digital convertor (ADC) and/or the digital-to-analog converter (DAC). In this work, an arbitrary waveform generator is set at 500 MSa/s, i.e., = 2ns to inject the ( ) • ( ) to the photonic reservoir. With a training data length of = 600 and = 14, the signal processing time is ~ 16.8 ( × × ). The full system computation time and the alternative approach in realizing the hardware-based RC for real-time implementation are discussed in the supplementary material.

Methods
The tumor respiratory motion data is collected from an external surrogate using a Varian TM system (RPM, Varian Medical Systems, Palo Alto, CA, USA). The sensor sampling rate is 30Hz. The radiator system latency in adjusting the multi-leaf collimator is caused by the combined software and hardware delay. In this work for the purpose of general study, we focus on short to medium range of equipment latency, i. e., the motion prediction look-ahead time of 66.6 ms, 166.6 ms and 333.3 ms.
Experimental setup. The real-time RC algorithm is implemented on a bench-top photonic TDR system constructed from off-the-shelf fiber optical components. A schematic representation of the TDR is shown in Fig. 8. The nonlinear node of the RC is realized by an optical amplitude modulator (EOSPACE AX-0MSS-20-PFA-SFA-LV). Adjusting the DC bias, the operation point of the modulator is set at 0.05 right to the peak of the sinusoidal transfer function, giving = 0.55 in equation (3). Besides several short fibers used to connect various optical components, there is no dedicated fiber used for extra time delay. The measured intrinsic time delay of the photonic TDR is 28 ns which comprises the true-time delay needed for the reservoir computer. A fiber-optic coupler (Newport F-CPL-1550-N-FA) and an electrical amplifier (Tektronix PSPL5865) are used in the TDR, the combining of which controls the feedback gain α.
It is worth noting that the electrical amplifier makes it easier to adjust the gain α, but it is not a must-have component and can be excluded in the photonic TDR system. The input gain is the coefficient of u(t) generated by an arbitrary waveform generator (AWG) (Keysight M8190), so it is only affected by the electrical portion of the delay loop (green path in Fig.8). The feedback gain is determined by both the optical path (red path in Fig.8) and the electrical path. The spike train method 27 presents visualization of the decaying memory of the RC and provides a simple approach to characterize how sensitive the reservoir dynamics in response to external stimuli. A decaying ratio is defined by the signal intensity between two consecutive pulses. In this work, a decay ratio of 0.87 is set [ Fig. 8 (c)]. The decaying ratio can also be treated as a hardware based hyperparameter and optimized for one task. The extracted feedback gain and values are = 0.87, = 2.29 based on measurement results. Details of and calculation are described in the supplementary material.