The approach is basically to use the restP attribute which controls the "default" position of the solid object, so that the mesh will try to keep it's initial shape, but that value could be animated and updated inside the sim using a sop solver. That restP attribute can be defined by vertex so that it can be transformed using the pivot of the tetrahedrons, so a "local transform" per each tetrahedron could be achieved, in this case contracting or expanding.
Then the rest is up to the FEM sim.
HIP file of one Slug with some comments inside:
A few videos i used for reference: