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.
However, as of version 11.1 specifying common CA rules has become a lot more convenient. The possibility to specify a CA rule via an association with allows for rather high-level classifications. In fact, Mathematica now knows about various neighbourhoods:
CellularAutomaton[<| "OuterTotalisticCode" -> 1018, "Neighborhood" -> "VonNeumann", |>, {#, 0} ][[1, 2 ;; -2, 2 ;; -2]] & And we don't even need to compute that rule code, because we can specify the rule directly via a set of growth cases:
CellularAutomaton[<| "GrowthCases" -> {2, 3, 4}, "Neighborhood" -> "VonNeumann", |>, {#, 0} ][[1, 2 ;; -2, 2 ;; -2]] & This says "when a dead cell has 2, 3 or 4 live neighbours, the cell comes alive", which is exactly what we're looking for.
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[<| "GrowthCases" -> {2, 3, 4}, "Neighborhood" -> "VonNeumann", |>, {#, 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:
To show the absolute age instead of the relative age, you can give ArrayPlot the option PlotRange -> {0, Length@b}:
