In optical communications, four-dimensional (4D) modulation formats encode information onto the quadrature components of two arbitrary orthogonal states of polarisation of the optical field. Many analytical models available in the optical communication literature allow, within a first-order perturbation framework, the computation of the average power of the nonlinear interference (NLI) accumulated in coherent fibre-optic transmission systems. However, all such models only operate under the assumption of transmitted polarisation-multiplexed two-dimensional (PM-2D) modulation formats, which only represent a limited subset of the possible dual-polarisation 4D (DP-4D) formats. Namely, only those where data transmitted on each polarisation channel are mutually independent and identically distributed. This paper presents a step-by-step mathematical derivation of the extension of existing NLI models to the class of arbitrary DP-4D modulation formats. In particular, the methodology adopted follows the one of the popular enhanced Gaussian noise model, albeit dropping most assumptions on the geometry and statistic of the transmitted 4D modulation format. The resulting expressions show that, whilst in the PM-2D case the NLI power depends only on different statistical high-order moments of each polarisation component, for a general DP-4D constellation, several other cross-polarisation correlations also need to be taken into account.