Membrane-mediated interaction between membrane-spanning peptides or protein segments plays an important role in their function and stability. Our rigorous "coupled harmonic oscillators" representation is extended to account for the complex boundary conditions permitting anisotropic relaxation of the membrane slope along the contours of the inclusions. Using this representation and applying a highly efficient finite-difference algorithm, we have analyzed the membrane-mediated interaction triggered by deformation of the hydrophobic tails of lipid molecules to match the lipophilic exterior of the inserted peptide. We establish that anisotropic relaxation crucially affects the interaction energy, leading to a short-range attraction between two inclusions, while conventional isotropic boundary conditions result in their strong repulsion. In a multi-inclusion cluster, this attraction is further enhanced and modified due to nonpairwise interactions. The results for dimyristoyl phosphatidylcholine and glyceryl monooleate membranes are compared, and the effects of the inclusion radius are considered. The possible role of slope relaxation in the reported stabilization of linked gramicidin channels and in proteins' functional cooperativity is outlined. (C) 2004 American Institute of Physics.