This paper deals with the oncolytic virotherapy model{ut=DuΔu−ξu∇⋅(u∇v)+μuu(1−u)−ρuuz,x∈Ω,t>0,vt=−(αuu+αww)v+μvv(1−v),x∈Ω,t>0,wt=DwΔw−ξw∇⋅(w∇v)−δww+ρwuz,x∈Ω,t>0,zt=DzΔz−δzz−ρzuz+βw,x∈Ω,t>0, which was initially proposed by Alzahrani-Eftimie-Trucu (2019) [1] to model the process of oncolytic viral therapy, which is a therapeutic approach for cancer treatment. Here, Ω⊆RN(N≥1) is a bounded smooth domain with zero-flux boundary conditions, where Du,Dw,Dz,δw,δz,μu,ρu are positive parameters and ξu,αu,ξw,ρw,ρz,μv are nonnegative constants. The main results assert the global boundedness of solutions to an associated spatially N-dimensional initial-boundary value problems under suitably assumptions on the system parameters. We develop (refine) some integral estimate techniques and thereby improve previous results in [18,23,14]. To the best of our knowledge, there are the first results on boundedness of the system in higher dimension (N≥2).