  This articles first investigates boundary integral operators for the
three-dimensional isotropic linear elasticity of a biphasic model with
piecewise constant Lam\'e coefficients in the form of a bounded domain of
arbitrary shape surrounded by a background material. In the simple case of a
spherical inclusion, the vector spherical harmonics consist of eigenfunctions
of the single and double layer boundary operators and we provide their spectra.
Further, in the case of many spherical inclusions with isotropic materials,
each with its own set of Lam\'e parameters, we propose an integral equation and
a subsequent Galerkin discretization using the vector spherical harmonics and
apply the discretization to several numerical test cases.
