Motion by mean curvature for networks of surfaces arises in a variety of applications, such as the dynamics of foam and the evolution of microstructure in polycrystalline materials. It is steepest descent (gradient flow) for an energy: the sum of the areas of the surfaces constituting the network. During the evolution, surfaces may collide and junctions (where three or more surfaces meet) may merge and split off in myriad ways as the network coarsens in the process of decreasing its energy. The first idea that comes to mind for simulating this evolution -- parametrizing the surfaces and explicitly specifying rules for cutting and pasting when collisions occur -- gets hopelessly complicated. Instead, one looks for algorithms that generate the correct motion, including all the necessary topological changes, indirectly but automatically via just a couple of simple operations. A remarkably elegant such algorithm, known as threshold dynamics, was proposed by Merriman, Bence, and Osher in 1992. Extending this algorithm, while preserving its simplicity, to more general energies where each surface in the network is measured by a different, possibly anisotropic, notion of area requires new mathematical understanding of the original version, which then elucidates a systematic path to new algorithms.