We analyze a plurality of epidemiological models through the lens of physics-informed neural networks (PINNs) that enable us to identify multiple time-dependent parameters and to discover new data-driven fractional differential operators. In particular, we consider several variations of the classical susceptible-infectious-removed (SIR) model by introducing more compartments and delay in the dynamics described by integer-order, fractional-order, and time-delay models. We report the results for the spread of COVID-19 in New York City, Rhode Island and Michigan states, and Italy, by simultaneously inferring the unknown parameters and the unobserved dynamics. For integer-order and time-delay models, we fit the available data by identifying time-dependent parameters, which are represented by neural networks (NNs). In contrast, for fractional differential models, we fit the data by determining different time-dependent derivative orders for each compartment, which we represent by NNs. We investigate the identifiability of these unknown functions for different datasets, and quantify the uncertainty associated with NNs and with control measures in forecasting the pandemic.