One thing I found surprising about this is that there are no random variables used, and the starting configuration is symmetrical. Yet it still buckles into these very asymmetrical shapes.

Floating point error accumulation ?

@Mateusz - indeed, this must be the cause

Could it be that the asymmetry comes from the imbalance caused by the sequential execution (iterating point by point)?

@Daniel - Kangaroo does not iterate point-by-point, but does a simultaneous global solve for each iteration.

I think it is more that although the starting configuration has icosahedral symmetry (it is a mesh subdivided with Weaverbird), when this is represented with finite precision xyz coordinates it introduces tiny numerical imperfections, which trigger the symmetry breaking.

What's nice is that slight variations in the input parameters (pressure, volume reduction fraction, smoothing strength, etc.), can produce such richly varying geometries.

Here's some of the research that prompted my interest in these matters:

http://web.mit.edu/preis/www/research.html

Ok, good to know! But where this imperfections may arise if the accuracy is the same for all points-coordinates? They should not have the same error at every point, so would maintain the symmetry? Or am I missing something? I saw the lines of the mesh after the subdivision are of different sizes, it is not most probable that local differences in the lengths of the lines triggered the broken symmetry? You'd be able to get those wrinkles in a mesh with identical edges as a box mesh?

Daniel: as Mateusz said it looks like it comes from floating point error accumulation.

Also helps that physical simulations of this complexity are really chaotic, so small floating point errors snowball into larger changes. In situations where a more stable equilibrium is reachable this isn't a concern.

http://en.wikipedia.org/wiki/Double_pendulum

I understand that, but assuming that idea, I wonder... having a spherically symmetric mesh (thi is, rebuilding the mesh projecting points to a spherical surface, all points equally far to the center), where the irregularity is born? Why a differential error should appear on floating values if all points are applied the same force? Rhino has a limit on the floating accuracy but should apply to all points in the same way? Why the floating error accumulation is not constant at all points?

I think the error are in the length of the lines of input mesh. Someone could get wrinkles asymmetrical in a input mesh with edges with the same length? I tried it and could not.
This is a unimportant issue, but I'm too curious to not deepen xD
Adding irregularities in edges, the mesh should wrinkle asymmetrically, and it is.

It seems that even in a spherical surface the floating points' "noise" is not the same for every point on the surface:

Measuring the distance from the center of the sphere to 100 random points I got 4 different numbers...