Details on specifying random parameters in 3-level HLM with lme

Dears,

I am used to fit 3-level HLM with Stata or HLM software, but now I have to
do it with lme, from package nlme. Unfortunatelly, though, I am struggling
with basic things and I could use any tips or suggestions.

What I am trying to do is to estimate a 3-level linear model where some of
random effects are in the 2nd level, some are in the 3rd. The specification
equations and the R code I was thinking of are as follows:

*Y* = b0 + b1**X1* + b2**X2* + e

      b0 = a00 + a01**X3* + r0
      b1 = a10 + a11**X4* + r1
      b2 = a20 + r2

              a00 = c00
              a01 = c10 + u1
              a10 = c20 + u2
              a11 = c30
              a20 = c40

*Y* = a00 + a01**X3* + r0 + a10**X1* + a11**X4***X1* + r1**X1* + a20**X2* +
r2**X2* + e

*Y* = c00 +c10**X3* + u1**X3* + r0 + c20**X1* + u2**X1* + c30**X4***X1* + r1
**X1* + c40**X2* + r2**X2* + e

*Y* = c00 + c20**X1* + c40**X2* + c10**X3* + c30**X4***X1* + r0 + r1**X1* +
r2**X2* + u2**X1* + u1**X3* + e

So far, I could only wonder this code:
lme(Y ~ X1 + X2 + X3 + X1:X4, random= ~ X1 + X2 + X3 | level3 / level2,
data=mydataset)

But I feel this is somewhat wrong, as it estimates the random parameters of
X1, X2 and X3 both in level 2 and in level 3, while the correct would be to
put the random compenents of X1 and x2 in level 2 and of X1 and X3 in level
3. Which is something I don't know how to separate in the "random"
parameter of the lme function and I couldn't find out anywhere so far.