You may be interested in these this paper

http://onlinelibrary.wiley.com/doi/10.1111/j.1538-4632.2000.tb00425.x/abstract

A Bayesian Approach to Prediction Using the Gravity Model, with an ...

and if you have nor seen it here are MLwiN resources on the cross-classified model

http://www.bristol.ac.uk/cmm/software/mlwin/mlwin-resources.html#nonhier

Although not specifically about cross-classified models, this post is worth reading

When Can You Safely Ignore Multicollinearity?

http://statisticalhorizons.com/multicollinearity

For the sake of illustration let's say you have a 400 x 400 system ... and that 75% of your 19310 observations are along the diagonal of this matrix -- i.e.from the same origin to destination. So, roughly 4820 observations are spread over the remaining 159600 cells!

Your data are very sparse as you can see. I would worry about HOW you are going to measure intra-area distance/travel time.

Then, I would want to know how many of the off-diagonal cells are structurally zero (i.e. impossible). Can you mask the ones that have no realistic chance of flow?

I think a little time playing with these concepts would help. What would help even more would be more data or a zoning scheme where the interactions are not all intra-zonal.

Finally, given the data you have, you might want to examine how you characterize the zone as an ORIGIN and as a DESTINATION.... since many (75%) of the flows will be from I to I.

From my experience with log linear models in SPSS I think you will have difficulty calibrating this kind of model with that software. 

A couple of things

1) what is your outcome variable and how is it measured?

2) Morton's post is useful but the nature of the problem changes according to weither you are using the additive cross classification or the multiplicative version; in the latter you have 424 BY 401  but in the former you would have just(! )  424 random effects PLUS 401 random effects which may be possible - that is there is an unexplained bit due to the origin place and due to the destination place but not due to their interaction. 

3) on the multicollinearity issue - I would also try and use a standard non-multilevel regression model with the variables in and look at the VIFs to see how bad the issue is. You could also create a use a variable that assesses the difference between the origin and destination values of the predictor and even say the origin value and the difference in the origin and destination values. I have never done this so I am just giving out ideas! 

Let us know how you get on.

Thank you both very much for the valuable advice, and apologies for the late response.

About my approach so far: my outcome variable is a binary indicator of mode choice (1=car, 0=other), to which I fit a logit model. I'm using an additive cross classification without random coefficients, so just the random effects by origin and destination area. MCMC shows acceptable mixing and effective sample size when I specify a large enough number of draws (I use 500,000 with thinning of 10). Fitting a variance components model suggests that about 95% of variance is attributable to the individual level, 3.3% to the origin area and 1.6% to the destination area (I understand from the CMM LEMMA course that using logit these VPC values are approximate rather than exact). 

I'm working on a method to move the analysis from the municipal scale to zipcode scale, which should more or less eliminate the issue of overlap between origin and destination areas. Or perhaps I should combine zipcode and municipality into a more complex cross-classified model, if this is still computationally feasible.