Discrete-continuous hybrid models are a popular means for describing elastic membrane-mediated particle interactions in and on lipid bilayers. Here, the continuous part is usually given by an approximation of the lipid membrane by an infinitely thin and sufficiently smooth hypersurface, whose elastic energy is determined by a Canham-Helfrich type functional. The discrete component results from modeling non-membrane particles as rigid discrete entities, which, depending on their configuration, induce local constraints on the membrane along the membrane-particle interfaces. In this context, the interaction potential describes the optimal elastic energy of such hybrid systems with a fixed particle configuration. Correspondingly, the energy minimization principle yields that stationary particle configurations are given by the local minima of the interaction potential. The main goal of this work is the proof of differentiability of the interaction potential for a selected class of models. This is accomplished using a variational approach that is already established in the literature in order to develop and apply robust numerical optimization methods for computing stationary particle configurations. Correspondingly, an additional focus is the derivation of a numerically accessible representation of the gradient, including its discretization and relevant numerical analysis. The proof of differentiability is brought forward by an application of the implicit function theorem. The basis for this is so-called boundary preserving domain transformations, which are induced by suitable families of vector fields and which locally admit the reformulation of the minimization problem that is implicitly defined by the interaction potential with respect to a fixed particle configuration. This subsequently enables the representation of the gradient as a volume integral using matrix analysis methods. The discretization of the partial differential equations for describing optimal membrane shapes is done via finite element methods. For particle methods with so-called curve restrictions a fictitious domain stabilized Nitsche method is developed, and for models with point value restrictions a conforming Galerking discretization is made possible by local QR transformations of the nodal finite element basis. For both cases suitable a priori error estimates are proven, and in addition also error estimates for the volume representation of the gradient are shown within that context. These developed methods open up the domain of efficient simulation of macro structures by isotropic and anisotropic particles, which is illustrated with the aid of various example applications and by means of perturbed gradient methods.