Hamiltonians of the superconducting qubits of transmon type involve nonzero ZZ interaction terms due to their finite and small anharmonicities. These terms might lead to the unwanted accumulation of spurious phases during the execution of the two-qubit gates. Exact calculation of the ZZ interaction rates requires the full diagonalization of the circuit Hamiltonians, which very quickly becomes computationally demanding as the number of the modes in the coupler circuit increases. Here we propose a direct analytical method for the accurate estimation of the ZZ interaction rates between low-anharmonicity qubits in the dispersive limit of the multimode circuit quantum electrodynamics (circuit QED). We observe very good agreement between the predictions of our method and the measurement data collected from the multiqubit devices.