Our focus is to study constellations of disjoint disks in the hyperbolic space, i.e., the unit disk equipped with the hyperbolic metric. Each constellation corresponds to a set E which is the union of m > 2 disks with hyperbolic radii r j > 0, j = 1, . . . , m. The centers of the disks are not fixed, and hence individual disks of the constellation are allowed to move under the constraints that they do not overlap and their hyperbolic radii remain invariant. Our main objective is to find computational lower bounds for the conformal capacity of a given constellation. The capacity depends on the centers and radii in a very complicated way even in the simplest cases when m = 3 or m = 4. In the absence of analytic methods, our work is based on numerical simulations using two different numerical methods, the boundary integral equation method and the hp-FEM method, respectively. Our simulations combine capacity computation with minimization methods and produce extremal cases where the disks of the constellation are grouped next to each other. This resembles the behavior of animal colonies minimizing heat flow in arctic areas.