If it is at all an option to represent the grid as a 2D list instead of a list of infected coordinates, I would model this is a cellular automaton. What you've essentially got is an outer totalistic cellular automaton with a von Neumann neighbourhood. The rule in Game-of-Life notation is B234/S01234, i.e. a cell comes to life if it has two or more live neighbours and it always survives. Implementing simple CAs is quite straight-forward with Mathematica's CellularAutomaton, and I've written another answer here about how to figure out the rule number of the CA.
For your case, we're using the weights:
{{0, 2, 0}, {2, 1, 2}, {0, 2, 0}} And then the rule turns out to be 1018. So we can simulate a single step with the following function:
CellularAutomaton[ { 1018, {2, {{0, 2, 0}, {2, 1, 2}, {0, 2, 0}}}, {1, 1} }, {#, 0} ][[1, 2 ;; -2, 2 ;; -2]] & The indexing at the end is used to remove the background information returned by CellularAutomaton.
To simulate the infection to convergence, I'd recommend FixedPointList instead of NestWhileList. It simply applies a function over and over until the value stops changing, and then gives you all the intermediate values.
Module[{a, b, d = 25}, a = RandomChoice[{0, 0, 0, 0, 0, 0, 0, 1}, {d, d}]; b = Most@ FixedPointList[ CellularAutomaton[{1018, {2, {{0, 2, 0}, {2, 1, 2}, {0, 2, 0}}}, {1, 1}}, {#, 0}][[1, 2 ;; -2, 2 ;; -2]] &, a]; ListAnimate[ArrayPlot /@ b] ] 
Adding some information about the history is as easy as calling Accumulate on the list of grids before handing them to ArrayPlot, which now colours each cell by its relative age:
