We present a lattice computation of pseudoscalar and vector heavy-light meson masses for heavy-quark masses ranging from the physical charm mass up to $\simeq 4$ times the physical b-quark mass. We employ the gauge configurations generated by the European Twisted Mass Collaboration (ETMC) with $N_f = 2+1+1$ dynamical quarks at three values of the lattice spacing ($a \simeq 0.062 - 0.089$ fm) with pion masses in the range $M_\pi \simeq 210 - 450$ MeV. Heavy-quark masses are simulated directly on the lattice up to $\simeq 3$ times the physical charm mass. The interpolation up to the physical $b$-quark mass is performed using the {ETMC ratio} method, based on ratios of (spin-averaged) meson masses and adopting the kinetic mass scheme. The extrapolation to the physical limit yields $m_b(\mbox{kin}, 1~\mbox{GeV}) = 4.61 (20)$ GeV, corresponding to $m_b(\overline{MS}, m_b) = 4.26 (18)$ GeV. Then the {\it ratio} method is applied above the physical $b$-quark mass to provide heavy-light meson masses towards the static point. The lattice data are analyzed in terms of the Heavy Quark Effective Theory (HQET) and the matrix elements of dimension-4 and dimension-5 operators are extracted with a good precision, namely: $\Lambda = 0.552 (26)$ GeV, $\mu_\pi^2 = 0.321 (32)$ GeV$^2$ and $\mu_G^2(m_b) = 0.253 (25)$ GeV$^2$. The data allows also to estimate the size of the matrix elements of dimension-6 operators. The extracted HQET parameters play a crucial role in the analysis of the inclusive semileptonic $B$-meson decays relevant for the determination of the CKM entry $V_{cb}$, this application is briefly discussed.