The superconductor-insulator transition (SIT) in thin-film disordered superconductors is a hallmark example of a quantum phase transition. Despite being observed more than 30 years ago, its nature is still under a vivid debate. One intriguing observations concerns the insulating side of the transition, which exhibits some unusual properties. Among them is its current-voltage relation ($I-V$ curve), which includes (1) a conductance that changes abruptly by several orders of magnitude with increasing voltage, (2) hysteretic behavior, and (3) multiple (sometimes more than a hundred) smaller current jumps near the transition. Some models have been suggested before, but no model has been successful in accounting for the observed behavior in full. One commonly used approach is to model the disordered sample as a two-dimensional array of conducting islands, where charge carriers tunnel from one island to its neighbors. In those models, fast relaxation is assumed, and the system is treated as always being in electrostatic and thermal equilibrium. Those models are successful in explaining some measurement results, including the phase transition itself, but fail to reproduce hysteresis in their predicted $I-V$ curves. Here, we suggest incorporating finite relaxation time into an array model. We show that, in the slow relaxation limit, our model can reproduce hysteresis and multiple jumps in the $I-V$ curve. Based on our results we argue that a similar behavior should be also observed in two-dimensional normal (non-superconducting) arrays. This claim is supported by past observations. We analyse the role of different parameters in our model, determine the range of relevant time scale in the problem, and compare our results with selected measurements.