Real-time implementation of cascaded all-pass filters from given transfer function

I am working on real-time implementation of spring reverb based on scientific paper called Parametric Spring Reverberation Effect by Välimäki, Vesa; Parker, Julian; Abel, Jonathan S.

One block of this effect is a cascade of $100$ stretched all-pass filters. Single stretched all-pass filter is defined by: $$A_{\rm low}(z) = \displaystyle\frac{a_1 + A_{\rm fd}(z)z^{-K_1}}{1 + a_1A_{\rm fd}(z)z^{-K_1}}=\displaystyle \frac{a_1 + \displaystyle\frac{a_2+z^{-1}}{1+a_2z^{-1}}z^{-K_1}}{1 + a_1\displaystyle\frac{a_2+z^{-1}}{1+a_2z^{-1}}z^{-K_1}}$$ I calculated the difference equation of this filter:

$$y(n)=a_1x(n)+a_1a_2x(n-1)+a_2x(n-K_1)+x(n-1-K_1)-a_2y(n-1)-a_1a_2y(n-K_1)-a_1y(n-1-K_1)$$

and used it for implementing the cascade of $M = 100$ stretched all-pass filters in MATLAB.

for i = 1:M
for n = 1+K1+1:(length(x))
y(n) = a1 * x(n) + a1*a2 * x(n-1) + a2 * x(n-K1) + ...
x(n-1-K1) - a2 * y(n-1) - a1*a2 * y(n-K1) - a1 * y(n-1-K1);
end
x = y;
end

However, the problem seems to be that I need to add another filters to the structure which are not cascaded (as you can see in the following block diagram), therefore they cannot be "repeated" like the all-pass filter. How can I calculate the output of my all-pass filter cascade sample-by-sample?

I can imagine the multiplication of two or three transfer functions in order to obtain the transfer function of the whole cascade, but multiplying $100$ transfer functions seems to be a crazy and wrong idea. There is probably a simple solution, but I just couldn't figure anything out. That's why I am asking you for help.