the system response by means of a circuit model, which includes an additional signal path introduced by the electromagnetic environment. Furthermore, we present two methods allowing us to experimentally determine the nonlinearity. First, we fit the measured frequency and flux dependence of the transmission data to simulations based on the equivalent circuit model. Second, we fit the power dependence of the transmission data to a model that is predicted by the nonlinear equation of motion describing the system. Our results show that we are able to tune the nonlinearity of the resonators by almost two orders of magnitude via an external coil and two on-chip antennas. The studied system represents the basic building block for larger systems, allowing for quantum simulations of bosonic many-body systems with a larger number of lattice sites.