Let $n(x,y)$ the number of agents in each location $(x,y)$ on the domain $S$.
The heterogeneous dynamics is described by a modulated Gaussian diffusion process.
Every time step $\textrm{d}t$, an agent moves from location $(x_0,y_0)$ to location $(x,y)$ with probability $\tilde{\mathrm{P}}(x,y|x_0,y_0)$. We refer to this as the transition kernel.
In the homogeneous case, the transition kernel $\mathrm{P}(x,y|x_0,y_0)$ is given by a 2D Gaussian \begin{equation} \mathrm{P}(x,y|x_0,y_0) = \frac{1}{2\pi\sigma^2}\exp\left(-\frac{(x-x_0)^2+(y-y_0)^2}{2\sigma^2}\right), \qquad (1) \end{equation} where $\sigma$ is the typical distance an agent moves in one timestep $\textrm{d}t$. That is, $\sigma=\sqrt{2 D \textrm{d}t}$ with $D$ the diffusion constant of the agent.
We modulate the homogeneous transition kernel by an affection map $a(x,y)$, representing the relative probability for an agent to move into (or stay in) position $(x,y)$. The modulated kernel then becomes \begin{equation} \tilde{\mathrm{P}}(x,y|x_0,y_0) = \frac{a(x,y) \mathrm{P}(x,y|x_0,y_0)}{\int_S a(x,y) \mathrm{P}(x,y|x_0,y_0) \textrm{d}x \textrm{d}y}. \qquad (2) \end{equation} Note that we renormalize this distribution, making the units of $a(x,y)$ arbitrary.
We now face the problem how to set the affection map $a(x,y)$ in order to achieve a certain steady state solution $\hat{n}(x,y)$.
Let $\textrm{d}n(x,y)$ the net number of agents moving into location $(x,y)$ during a certain timestep. Then, its expectation value is given by \begin{equation} \langle\textrm{d}n(x,y)\rangle = \underbrace{-n(x,y)}_{\textrm{outflux}} + \underbrace{\int_S n(x_0,y_0) \tilde{\mathrm{P}}(x,y|x_0,y_0) \textrm{d}x_0 \textrm{d}y_0}_{\textrm{influx}}. \qquad (3) \end{equation}
The steady state distribution $\hat{n}(x,y)$ solves the equation $\langle\textrm{d}n(x,y)\rangle = 0$. Hence, the affection map $a(x,y)$ that results in a given steady state solution $\hat{n}(x,y)$ can be obtained by solving \begin{equation} -\hat{n}(x,y) + \int_S \hat{n}(x_0,y_0) \tilde{\mathrm{P}}(x,y|x_0,y_0) \textrm{d}x_0 \textrm{d}y_0 = 0, \qquad (4) \end{equation} with eqs (1,2) for $a(x,y)$.
I have solved equation (4) numerically for a discretized case, but I wonder if it is possible to come up with an analytical solution of equation (4) for $a(x,y)$ through functional calculus. An additional equation would be required to set the scale of $a(x,y)$, so say $a(0,0)=1$.