The solution of this ODE defined over a compact interval $\left[0,T\right]$ for each initial condition can be formally cast as the limit (as the number of individuals $N$ scales to infinite) of a stochastic process modelling the peer-to-peer spread of a virus with Poisson clocks.
Remark. This answer is just intended to provide some references on the matter and break down a bit the overall standard path (along these references) that show how these non-stochastic ODEs emerge formally from the stochastic interacting dynamics among $N$ individuals (or nodes in a complete graph) in the thermodynamic limit, i.e., as $N\longrightarrow \infty$.
In words. i) characterize the microscopic state of the system (i.e., the state of each of the $N$ individuals as infected, susceptible or removed over time); ii) via peer-to-peer interacting rules, ascertain the dynamics, i.e., how the microscopic state will evolve over time according to neighbouring interactions (one may resort to an interacting particle system [1], where the interactions among individuals are driven by Poisson clocks); iii) Characterize the macroscopic processes, e.g., the fraction of infected individuals or the fraction of susceptible individuals over time (these processes build on the microscopic one); iv) take the limit as $N\rightarrow \infty$ and under certain conditions the macroscopic processes converge to the solution of the SIR or SIS ODE (e.g., following Kurtz approach [2-3], the macro processes are governed by the solution to the ODE up to a martingale term that vanishes as $N$ grows to infinite, yielding the concentration).
Just to depict an overall thermodynamic limit program in a nutshell, you have the following possibility.
Microscopic process. Consider $N$ individuals. Let $X^N_i(t)\in\left\{-1,0,1\right\}$ to be a process representing the state of an individual $i$ at time $t\geq 0$: i) $X^N_i(t)=0$ means that individual $i$ is susceptible; ii) $X^N_i(t)=1$ means infected; and iii) $X^N_i(t)=-1$ means removed at time $t$.
Microscopic dynamics [1]. Once infected, each individual $i$ turns on two clocks: a) Infection clock -- after a random time $T_I\sim {\sf Exp}(1/a)$, it randomly chooses another individual and infects him -- if the chosen individual is already infected, then nothing happens; b) Removed clock-- after a random time $T_r\sim {\sf Exp}(1/b)$, it transitions to the state removed (and the clocks are turned off).
Macroscopic process. Define $I^N(t):= \sum_{i=1}^N \mathbf{1}_{\left\{X^N_i(t)=1\right\}}/N$; $S^N(t):= \sum_{i=1}^N \mathbf{1}_{\left\{X^N_i(t)=0\right\}}/N$ and $R^N(t):= \sum_{i=1}^N \mathbf{1}_{\left\{X^N_i(t)=-1\right\}}/N$ as the fraction of infected, susceptible and removed individuals at time $t$ among $N$ individuals in the population.
Macroscopic dynamics. This characterization yields -- the macro-processes will be characterized by the integral of the Poisson processes governing the microscopic interactions and via compensating these Poisson processes, you obtain --
$$S^{N}(t) = S^{N}(0) + M_S^N(t) - a\int_{0}^t S^N(s-)I^N(s-)ds$$
$$I^{N}(t) = I^{N}(0) + M_I^N(t) + a\int_{0}^t S^N(s-)I^N(s-)ds-b\int_{0}^t I^N(s-)ds,$$
where $M_S^N(t)$ and $M_I^N(t)$ are martingales that converge in probability to zero (w.r.t. the uniform norm on the space of Skorokhod paths defined over the compact interval $\left[0,T\right]$, $\mathcal{D}\left[0,T\right]$). This will further imply that the processes $\left(I^N(t)\right)_{t\in\left[0,T\right]}$ and $\left(S^N(t)\right)_{t\in\left[0,T\right]}$ converge weakly (w.r.t. the uniform norm on $\mathcal{D}\left[0,T\right]$) to the solution of the ODE with initial condition $I(0)$ and $S(0)$ if $I^N(0)$ and $S^N(0)$ converge in distribution to $I(0)$ and $S(0)$.
More on the convergence. That the martingale converges to zero is granted by Doob's inequality. This will imply that the stochastic macro-processes $\left(I^N(t)\right)_{t\in \left[0,T\right]}$ and $\left(S^N(t)\right)_{t\in \left[0,T\right]}$ are tight (a useful result in this part is the stochastic version of Arzelà–Ascoli theorem, e.g., in [4]). This implies that $\left(I^N(t)\right)_{t\in \left[0,T\right]}$ and $\left(S^N(t)\right)_{t\in \left[0,T\right]}$ admit convergent subsequences whose limit is necessarily the solution of the ODE. Since these ODE's have unique solution (vector field is Lipschitz and they live in a compact invariant set $\left[0,1\right]^2$) the sequence converges.
Remark. If the graph of interactions is not complete, then $\left(I^N(t)\right)_{t\in \left[0,T\right]}$ and $\left(S^N(t)\right)_{t\in \left[0,T\right]}$ are not Markov and the problem becomes quite intricate. I believe that there is now a growing body of literature devoted to establishing this limit over nontrivial non-complete graphs.
[1] T. M. Liggett, Interacting particle systems, Springer, New York, 1985.
[2] T. G. Kurtz, Solutions of ordinary differential equations as limits of pure jump markov processes, Journal of Applied Probability, 7 (1970), pp. 49–58.
[3] S. N. Ethier and T. G. Kurtz, Markov Processes: Characterization and Convergence, ser. Probability and Statistics. New York, NY: John Wiley & Sons, Inc., 1986.
[4] J. Jacod and A. Shiryaev, Limit Theorems for Stochastic Processes, 2nd ed., ser. Grundlehren der mathematischen
Wissenschaften. New York, NY: Springer-Verlag, 2003, vol. 288.
