What happens when this integrator circuit reaches saturation?
What is the formula for calculating the time taken for V(R3) to drop to zero? Is it just 5RC (which is 5ms) and tallies with the simulation or do I have to take into account the fact that there's charge already accumulated on C1 plates due negative voltage at Vout?


It's an exponential decay function, so it never really reaches zero. You just have to decide how close is "close enough". 5&times;RC is the time that it takes to be within 1% of its final value. 3&times;RC gets you within 5%. The initial charge has no appreciable effect on this.

The general formula is this: The time it takes for the exponential decay curve to be within 1/N (we're talking about a fraction of the difference between the initial and final values) of the final value is

$$t = \ln(N)\cdot\tau$$

where tau (\$\tau\$) is the "time constant", equal to RC in this case. So, ln(20) = 3 and ln(100) = 4.6, which we generally round up to 5.


  Additionally, what happens when Vin changes to -5V (at 21ms). What's the formula for calculating the time taken for V- to reach zero? 


This is just the first part of the exponential decay, where the value reaches the point halfway between its initial and final values. As anyone who has used a 555 timer knows, this is 0.693&times;RC. (Because ln(2) = 0.693147)