I'm not sure if this is what you're looking for, but you can always describe the evolution of a many-boson system as the evolution of a set of qudits in the standard "kets and unitaries" formalism, at least provided the evolution is linear (as in, the total number of bosons is preserved). Probably also more generally, but I'm not sure.
Simple example
Consider a simple example. Two bosons in two modes, evolving with a unitary $U$. Say the initial state is $a_1^\dagger a_2^\dagger |0\rangle\equiv |1_1 1_2\rangle$, meaning one boson in each input mode. If the bosons evolve through the (single-particle) unitary $U$, then the overall output state is obtained evolving each creation operator via $U$, meaning applying the mapping
$$a_i^\dagger \to \sum_j u_{ji} a_j^\dagger,$$
and thus reads
$$|\operatorname{output}\rangle= (u_{11} a_1^\dagger + u_{21} a_2^\dagger)(u_{12} a_1^\dagger + u_{22} a_2^\dagger)|0\rangle
= \Bigg[\sqrt2 u_{11} u_{12} \frac{(a_1^\dagger)^2}{\sqrt2} + (u_{12}u_{21}+u_{11}u_{22})a_1^\dagger a_2^\dagger + \sqrt2 u_{21} u_{22} \frac{(a_2^\dagger)^2}{\sqrt2} \Bigg] |0\rangle
= \sqrt2 u_{11} u_{12} |2_1 0_2\rangle + (u_{12}u_{21}+u_{11}u_{22})|1_1 1_2\rangle + \sqrt2 u_{21} u_{22} |0_1 2_2\rangle.$$
But there's another way to describe similar situations. If we were not using the second quantisation formalism (that is, no creation/annihilation operators), then we would describe the initial state of two bosons in the two input modes with something like $|1,2\rangle$, where now the first (second) number tells you the state of the first (second) boson. But the bosons are indistinguishable aside for their mode state, thus we cannot say things like "the first boson is in the first mode and the second in the second mode", we have to say things like "one boson is in the first mode and the other one in the second mode". We can do this by symmetrising the input state, with something like
$$\frac1{\sqrt2}(|1,2\rangle+|2,1\rangle).$$
Now each boson evolves independently with $U$, hence the output state reads
$$(U\otimes U)\frac1{\sqrt2}(|1,2\rangle+|2,1\rangle).$$
You might now see that this allows you to simulate the situation with a two-qubit system. Upon redefinition of the labels, this is a Bell state evolving through a local unitary $U\otimes U$. There's however a catch, that you'll also want to perform non-local measurements on the output. This is because in boson sampling you want the output probabilities associated to the various occupation numbers, and the three (in this simple example) corresponding basis states read in this notation:
$$|2_1 0_2\rangle \sim |1,1\rangle,
\qquad |0_1 2_2\rangle \sim |2,2\rangle,
\qquad |1_1 1_2\rangle \sim |1,2\rangle+|2,1\rangle.$$
You can verify explicitly that doing this produces exactly the same probability amplitudes the original description of the many-boson system did.
General statement
The above is fully general. Given an $n$-boson, $m$-mode state, you can describe it as a symmetrised $n$-qudit state with each qudit having dimension $m$. Linear evolution is then described by some $U\otimes\cdots\otimes U$, and the outcome probabilities for the many-boson system can be recovered projecting on the corresponding $n$-qudit states. This sort of thing you can do in a circuit, and I suppose if you can simulate $n$-qudit states with $nm$-qubit ones, you can build your circuit with qubits.
Not that this will be easy, mind you. I'd expect the hardness of boson-sampling will translate into the hardness of creating the highly entangled states necessary to play this little game.
