In this paper, a comprehensive integral equation formulation of plasmonic transmission lines is presented for the first time. Such lines are made up of a number of metallic strips with arbitrary shapes and dimensions immersed within a stack of planar dielectric or metallic layers. These lines support a number of propagating modes. Each mode has its own phase constant, attenuation constant, and field distribution. The presented integral equation formulation is solved using the Method of Moments (MoM). It provides all the propagation characteristics of the modes. The new formulation is applied to a number of plasmonic transmission lines, such as: single rectangular strip, horizontally coupled strips, vertically coupled strips, triangular strip, and circular strip. The numerical study is performed in the frequency (wavelength) range of 150-450 THz (0.66-2.0 μm). The results of the proposed technique are compared with those obtained using Lumerical mode solution, and CST. Very good agreement has been observed. The main advantage of the MoM is its intrinsic speed for this type of problem compared to general purpose solvers.