The careful analysis of strongly gravitationally lensed radio and optical images of distant galaxies can in principle reveal DM (sub-)structures with masses several orders of magnitude below the mass of dwarf spheroidal galaxies. However, analyzing these images is a complex task, given the large uncertainties in the source and the lens. Here, we leverage and combine three important computer science developments to approach this challenge from a new perspective. (a) Convolutional deep neural networks, which show extraordinary performance in recognizing and predicting complex, abstract correlation structures in images. (b) Automatic differentiation, which forms the technological backbone for training deep neural networks and increasingly permeates 'traditional' physics simulations, thus enabling the application of powerful gradient-based parameter inference techniques. (c) Deep probabilistic programming languages, which not only allow the specification of probabilistic programs and automatize the parameter inference step, but also the direct integration of deep neural networks as model components. In the current work, we demonstrate that it is possible to combine a deconvolutional deep neural network trained on galaxy images as source model with a fully-differentiable and exact implementation of the gravitational lensing physics in a single probabilistic model. This does away with hyperparameter tuning for the source model, enables the simultaneous optimization of nearly one hundred source and lens parameters with gradient-based methods, and allows the use of efficient gradient-based Hamiltonian Monte Carlo posterior sampling techniques. We consider this work as one of the first steps in establishing differentiable probabilistic programming techniques in the particle astrophysics community, which have the potential to significantly accelerate and improve many complex data analysis tasks.