Start from the 1st law of thermodynamics $dE=\delta Q+\delta W$. Forget the work; this tells you that if heat flows into the system, its internal energy will change accordingly.
But this is a global view: it doesn't tell anything about where the heat flows through.
If you refine the description, you will define a local heat flow everywhere on the boundary of the system. This will of course be a vector : the heat can flow into the system, or out of the system, or in any direction really - if the heat comes from a temperature difference, it will flow towards lower temperatures. Let's call this local heat flow $\mathbf q$. The projection of this flow onto an unit vector $n$ is $\mathbf q \cdot \mathbf n$; it represents be the heat flow in the direction pointed to by n (note that if n is perpendicular to the heat flow, no heat enters the system as expected).
Now for an infinitely small surface element of area $da$ on the boundary of the system, the heat entering the system per unit time will be $- \mathbf q \cdot \mathbf n da$ where $\mathbf n$ is the outward unit vector normal to the surface. Integrate over the whole surface and you get the heat recieved by the system through the whole surface, per unit time : $$\frac{\delta Q}{dt} = - \int_A \mathbf q \cdot \mathbf n da$$
Then, neglecting work, the 1st law of thermodynamics gives you the wanted formula.
This is actually a pretty standard calculation for conserved quantities (that is, stuff that we don't expect to appear out of nowhere, such as mass, electric charge, energy, and many others). The variation of the total amount of stuff in your system is expected to be related to the flux of stuff crossing the system boundary, plus the stuff created inside the system. This leads to the equation $$\frac{dM}{dt}=-\int_A \mathbf Γ_m\cdot \mathbf n da + S_m$$ where $M$ is the total amount of stuff in the system, $\mathbf Γ_m$ is the flux of stuff, and $S_m$ is the amount of stuff created inside the system; $M$ is related to the density of stuff $m$ by $M=\int_V mdV$. This is equivalent to the local PDE $$\frac{\partial m}{\partial t}+\mathbf\nabla \mathbf Γ_m=s_m$$ familiar in fluid dynamics amongst many other fields.